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ABSTRACT 


This  thesis  explores  Benders  decomposition  for  solving  interdiction 
problems  on  electric  power  grids,  with  applications  to  analyzing  the  vulnerability 
of  such  grids  to  terrorist  attacks.  We  refine  and  extend  some  existing 
optimization  models  and  algorithms  and  demonstrate  the  value  of  our  techniques 
using  standard  reliability  test  networks  from  IEEE. 

Our  implementation  of  Benders  decomposition  optimally  solves  any 
problem  instance,  in  theory.  However,  run  times  increase  as  Benders’  cuts  are 
added  to  the  master  problem,  and  this  has  prompted  additional  research  to 
increase  the  decomposition’s  efficiency.  We  demonstrate  empirical  speed  ups 
by  dropping  slack  cuts,  solving  a  relaxed  master  problem  in  some  iterations,  and 
using  integer  but  not  necessarily  optimal  master-problem  solutions.  These  mixed 
strategies  drastically  reduce  computation  times.  For  example,  in  one  test  case, 
we  reduce  the  optimality  gap,  and  the  time  that  it  takes  to  achieve  this  gap,  from 
16%  in  75  hours  to  5%  in  16  minutes. 
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DISCLAIMER 


The  reader  is  cautioned  that  computer  programs  developed  in  this 
research  may  not  have  been  exercised  for  all  cases  of  interest.  While  every 
effort  has  been  made  to  ensure  that  the  programs  are  free  of  computational  and 
logic  errors,  they  cannot  be  considered  fully  validated.  Any  application  of  these 
programs  without  the  additional  verification  is  at  the  risk  of  the  planner. 
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EXECUTIVE  SUMMARY 


This  thesis  extends  existing  optimization  models  and  methods  for 
analyzing  the  vulnerability  of  electric  power  grids  to  terrorist  attacks. 

Electric  power  systems  are  critical  to  the  United  States’  economy  and 
security.  After  the  events  of  September  1 1 ,  2001 ,  and  the  United  States’  ensuing 
war  against  terrorism,  the  fear  of  retaliation  against  critical  infrastructures  has 
become  a  major  concern  for  security  analysts.  The  vulnerability  of  the  electric 
power  systems  to  physical  disruptions  heads  the  list  of  concerns,  because  the 
U.S.  transmission  grid  has  not  expanded  as  quickly  as  demand  has  over  the  last 
decade,  and  thus  the  system  has  become  “brittle.”  This  brittleness  (vulnerability) 
became  more  evident  on  14  August  2003  when  a  relatively  modest  number  of 
system  malfunctions  caused  a  blackout  of  the  northeastern  region  of  the  country, 
and  raises  this  question:  How  much  worse  might  the  blackout  have  been  had  it 
resulted  from  a  well-planned  terrorist  attack? 

This  thesis  first  considers  certain  modeling  issues  of  the  problem  of 
optimal  interdiction,  and  then  focuses  on  techniques  for  solving  the  proposed 
models.  This  should  help  us  identify  the  most  critical  components  of  the  U.S. 
power  grid,  and  ultimately  help  determine  effective  protective  measures  to 
improve  the  system’s  robustness. 

We  refine  and  extend  the  existing  optimization  models  and  algorithms  of 
Salmeron  et  al.  that  identify  critical  system  components  (e.g.,  transmission  lines) 
by  creating  maximally  disruptive  attack  plans  for  terrorists,  who  are  assumed  to 
have  limited  offensive  resources.  Most  importantly,  this  thesis  exploits  bi-level 
structures  embedded  in  the  interdiction  models  to  enable  faster  solutions  through 
the  use  of  Benders  decomposition. 

We  first  validate  the  DC  power-flow  model  that  underlies  in  the  interdiction 
models.  By  manipulating  a  few  equations,  the  inherent  nonlinearities  in  these 
models  are  substituted  by  linear  forms,  converting  those  models  into  (linear) 

xvii 


mixed-integer  problems  (MIPs).  We  explore  Benders  decomposition  for  solving 
the  MIPS,  and  successfully  apply  it  to  two  standard  reliability  test  networks  from 
IEEE. 

Although  we  solve  the  IEEE  test  cases  to  optimality,  we  notice  that  the 
efficiency  of  the  Benders’  master  problem  deteriorates  as  (a)  the  number  of 
components  in  the  test  cases  increases,  and  (b)  the  number  of  Benders’  cuts  in 
the  master  problem  grows.  We  investigate  several  techniques  in  order  to 
improve  the  algorithm’s  speed.  First,  we  demonstrate  faster  convergence  by 
solving  the  mixed-integer  master  problem  exactly  every  iteration  only,  and 
solving  its  easier,  linear-programming  relaxation  otherwise.  Importantly,  the 
algorithm  thus  modified  always  maintains  a  valid  upper  (optimistic)  bound,  just  as 
the  original  does.  Since  relaxed  (non-integer)  master-problem  solutions  cannot 
be  used  by  the  subproblems  to  compute  lower  (pessimistic)  bounds,  we  exploit 
sub-optimal  integer  solutions  to  the  master  problem  every  iteration  to  improve 
the  lower  bound.  In  one  of  our  two  test  networks,  this  combined  techniques 
reduce  computational  time  by  one  third.  We  also  show  that  limiting  the  number 
of  cuts  in  the  master  problem  can  speed  convergence,  and  study  several 
strategies  for  keeping  or  deleting  cuts.  Results  are  mixed,  but  improvements  can 
be  dramatic.  For  example,  in  the  second  of  our  test  networks,  by  limiting  the 
number  of  cuts  to  500  and  solving  a  relaxed  master  problem  nine  of  every  ten 
iterations,  we  achieve  an  optimality  gap  under  5%  in  16  minutes;  this  compares 
to  a  16%  gap  in  75  hours  using  the  original  decomposition  algorithm 


GLOSSARY 


The  following  briefly  defines  the  general  concepts  related  to  electric  power 
networks  and  interdiction.  Some  of  these  definitions  are  from  (or  strongly  based 
on)  the  glossaries  provided  by  Energy  Information  Administration  (EIA  2003), 
Elec-Saver  (2003)  and  SIEMENS(2004). 

AC  Power:  Power  associated  with  alternating  current  circuits. 

Admittance:  Property  that  allows  the  flow  of  electrical  current  through 
reactive  circuit  elements  under  the  action  of  a  potential  difference.  Admittance  is 
the  reciprocal  of  impedance.  Admittance  equations  establish  the  relationship 
between  current,  impedance  and  voltage. 

Alternating  Current:  Currenf  that  periodically  reverses  direction. 

Bus  (or  Busbar):  A  heavy,  rigid  electrical  conductor  that  makes  a  common 
connection  between  several  electrical  circuits. 

Capacitance:  The  property  of  a  circuit  that  allows  it  to  store  an  electrical 
charge. 

Case:  A  set  of  data  to  be  analyzed,  along  with  the  results  of  the  analysis. 
Data  consists  of  power  grid  components  {lines,  buses,  generators,  substations, 
etc),  physical  data  (i.e.,  such  as  line  impedances,  generating  capacities,  etc.), 
non-physical  data  (e.g.,  interdiction  resource,  optimization  parameters,  etc.). 
Results  include  the  interdicted  components,  associated  generation  and  power 
flows,  and  load  shedding  for  every  period,  among  others. 

Costumer  Sector:  A  type  of  load  with  specific  requirements  (e.g.,  amount 
of  power  demand  and  cost  for  failing  to  provide  it). 

Current:  The  flow  of  electrons  in  a  circuit.  Current  is  measured  in 
amperes. 

DC  Power:  Power  associated  with  direct  current  circuits. 
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Direct  Current:  Current  that  flows  in  one  direction 


Disruption:  The  cost  of  power  shed  (in  dollars  per  hour)  or  the  cost  of 
energy  shed  (in  dollars),  as  a  consequence  of  interdiction. 

Energy:  The  result  of  integrating  power  over  time. 

Energy  Shed:  Amount  of  energy  that  cannot  be  supplied  to  the  load  (one 
or  several  customer  sectors)  over  the  course  of  a  given  period. 

Generation  (of  electricity):  The  production  of  electric  energy  through  the 
transformation  of  other  forms  of  energy.  The  amount  of  electric  energy 
produced. 

Generating  Unit  (or  Generator):  Any  combination  of  physically  connected 
generator(s),  reactor(s),  boiler(s),  combustion  turbine(s),  or  other  prime  mover(s) 
operated  together  to  produce  electric  power. 

Impedance:  The  total  opposition  to  alternating  current.  Impedance  is  the 
vector  sum  of  resistance  and  reactance.  The  unit  for  impedance  is  the  ohm. 

Inductance:  The  property  of  an  electrical  circuit  that  causes  it  to  oppose 
changes  in  current. 

Interdiction  Plan:  Specific  subset  of  the  electric  system  equipment  that 
might  be  interdicted  by  terrorists.  Optimal  or  near-optimal  interdiction  plans  are 
identified  in  for  given  interdiction-resource  scenarios. 

Interdiction  Resource:  A  numerical  value  associated  with  a  mathematical 
expression  that  represents  the  capacity  of  terrorists  to  carry  out  attacks.  For 
example  if  such  an  expression  has  the  form:  “(3  x  total  number  of  attacks  to 
buses)  +  (1  X  total  number  of  attacks  to  lines)  <  5,”  then  “5”  is  the  interdiction 
resource. 

Line:  See  “Transmission  Line.” 

Load:  Demand  for  electric  power,  measured  in  watts,  at  a  specific  point  in 

time. 
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One-line  Diagram:  Schematic  drawing  of  an  electrical  power  system  that 
uses  graphical  symbols  to  represent  electrical  equipment  such  as  buses, 
generators,  loads,  transmission  lines  and  transformers.  It  may  incorporate 
numerical  values  for  the  system,  such  as  line  power  flows,  generating  unit 
outputs,  bus  voltages,  etc. 

Phase  angle:  The  angle  by  which  the  sine  curve  of  the  voltage  in  a  circuit 
element  (or  a  combination  of  elements)  leads  or  lags  the  sine  curve  of  the  current 
in  that  circuit  element(s). 

Period  (of  restoration):  Each  of  the  stages  that  an  electric  power  network 
undergoes  following  an  attack,  as  interdicted  components  are  repaired  or 
replaced  over  time. 

Power  (transmission):  The  transfer  rate  for  energy.  Unit  of  measurement  is 
watts.  Power  is  sometimes  used  loosely  to  refer  to  electricity. 

Power  Shed:  Amount  of  power  that  cannot  be  supplied  to  a  load  or  loads 
(one  or  several  customer  sectors)  at  a  specific  point  in  time. 

Reactance:  Opposition  to  the  flow  of  alternating  current.  Measured  in 

ohms. 

Reactive  Power:  Power  associated  with  inductance  or  capacitance. 

Resistance:  The  propensity  of  a  circuit  to  oppose  current  flow.  It  is 
measured  in  ohms. 

Scenario:  A  particular  value  of  the  interdiction  resource.  In  our 
formulation,  we  analyze  the  worst-possible  Interdiction  Plans  for  a  given 
scenario. 

Slack  Bus:  (Also  called  swing  bus.)  A  bus  whose  (initially  unspecified) 
power  generation  must  be  determined  in  order  to  match  supply  and  demand  in 
the  network. 

Substation:  Facility  with  equipment  that  switches,  changes,  or  regulates 
electric  voltage  and  current. 
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Transformer:  A  static  electrical  device  which,  by  electromagnetic 

induction,  regenerates  AC  power  from  one  circuit  into  another  and/or  changes 
the  voltage  of  alternating  current. 

Transmission;  The  movement  or  transfer  of  electric  energy  over  an 
interconnected  group  of  lines  and  associated  equipment  between  points  of 
supply  and  points  at  which  it  is  transformed  for  delivery  to  consumers,  or  is 
delivered  to  other  electric  systems.  Transmission  is  considered  to  end  when  the 
energy  is  transformed  for  distribution  to  the  consumer. 

Transmission  Line  (High-voltage):  The  high-voltage  (>  69  kV  in  the  U.S.) 
conductors  used  to  carry  electrical  energy  from  one  location  to  another. 

Transmission  System  (also  referred  to  as  Electric  Power  Grid  and  Electric 
Network  in  this  thesis):  An  interconnected  group  of  electric  transmission  lines 
and  associated  equipment  for  moving  or  transferring  electric  energy  in  bulk 
between  points  of  supply  and  points  at  which  it  is  transformed  for  delivery  over  a 
lower-voltage  distribution  system  to  consumers,  or  is  delivered  to  other  electric 
systems. 


I.  INTRODUCTION 


This  research  continues  the  development  of  mathematical  models  and 
optimization  methods  for  planning  electrical  power  grids  that  are  robust  to 
terrorist  attacks.  We  focus  on  solving  an  interdiction  model,  which  is  represented 
as  a  mixed-integer  program  (MIP),  through  the  use  of  Benders  decomposition. 


A.  THE  ELECTRIC  POWER  GRID  IN  THE  U.S. 

The  electric  power  grid  in  the  United  States  (U.S.)  consists  of  over  10,000 
generating  units  with  a  total  production  in  excess  of  760  giga-watts  (GW),  and 
over  700,000  miles  of  transmission  lines,  all  controlled  by  approximately  100 
control  centers.  (Actually,  only  modest  interconnection  capacity  exists  between 
three  main  sub-systems  in  this  grid,  namely  the  Eastern  interconnection,  the 
Western  interconnection  and  the  Texas  interconnection.  Therefore,  “the  U.S. 
power  grid”  may  be  viewed  as  three,  nearly  independent  systems.) 

In  the  past,  utilities  have  been  concerned  with  the  grid’s  vulnerabilities  to 
isolated  natural  disasters  and  unscheduled,  technical  outages  caused  by 
equipment  failures.  To  date,  protecting  the  electric  grid  against  multiple,  nearly 
simultaneous  failures  has  not  been  a  high  priority  for  utilities  because  of  the 
expense  involved  and  the  relatively  low  frequency  of  such  events.  However,  in 
the  current  era,  a  set  of  nearly  simultaneous  failures  might  be  the  precise 
objective  of  a  terrorist  attack.  “Low  frequency”  may  no  longer  be  so  slow. 

As  a  result  of  the  threat  that  terrorism  currently  poses  to  critical 
infrastructure  of  all  types  in  the  U.S.,  understanding  the  physical  vulnerability  of 
the  electric  grid  has  become  more  important  not  only  for  the  electric  power 
industry,  but  also  for  national  security  (NSTAC  1997).  In  fact,  a  National  Security 
Assessment  (1997)  states  that  man-made  physical  destruction  is  the  greatest 
threat  facing  the  electric  power  infrastructure. 

The  U.S.  electric  power  grid  is  plagued  by  a  number  of  inherent 
weaknesses.  Among  the  most  significant  of  these  weaknesses  are  the  traversals 
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of  sparsely  populated  areas,  which  readily  expose  the  system  to  malicious 
attacks,  and  the  decreasing  reserve  in  transmission  capacity. 

The  U.S.  transmission  network  is  built  with  some  redundancy  and  multiple 
protective  devices;  however,  current  reserve  levels  in  transmission  capacity  may 
be  insufficient  to  back  up  multiple  failures  in  some  key  components,  as 
exemplified  by  the  14  August  2003  power-system  failure  (which  cost  our  national 
economy  between  $7  and  $10  billion  dollars;  ICF  Consulting  2003).  Not 
surprisingly,  this  failure  brings  additional  speculation  about  terrorist  threats  to  the 
U.S.  power  grid  (Germain  2003),  and  raises  important  questions  about  the 
resilience  of  the  U.S.  power  grid  and  the  possibility  of  widespread  loss  of  power 
to  customers. 

The  August  14*^^  blackout  could  be  financially  insignificant  when  compared 
to  the  effects  of  a  coordinated  terrorist  attack  on  several  key  facilities  during  a 
period  of  peak  load.  ICF  Consulting  (2003),  which  analyzed  the  cost  of  the 
August  14*^^  blackout,  states  that  a  terror-induced  blackout  could  prove 
significantly  more  costly  and  have  debilitating  impacts  on  the  affected  region  as 
well  as  the  entire  country.  Additionally,  a  simulated  terrorist  attack  on  the  Pacific 
Northwest's  power  grid  was  recently  conducted  under  the  project  name  “Blue 
Cascades.”  An  analysis  of  that  simulated  attack  showed  that  a  real  attack,  if 
successful,  could  wreak  havoc  on  the  nation's  economy,  shutting  down  power 
and  productivity  in  a  domino  effect  that  would  last  for  weeks  (Allen  2003). 

The  discussion  above  motivates  the  development  of  optimization  models 
to  represent  the  problem  of  terrorists  attacking  a  power  grid.  By  studying  how  to 
attack  power  grids,  one  can  gain  insight  on  how  to  protect  them. 

This  research  follows  Salmeron  et  al.  (2003,  2004),  who  develop  bi-level 
mathematical  models  to  identify  maximally  disruptive  attacks  for  terrorists,  who 
are  assumed  to  have  full  access  to  power-system  data,  but  limited  offensive 
resources.  Salmeron  et  al.  initially  introduce  a  static  interdiction  model  and  then 
extend  it  to  capture  the  dynamics  of  system  operation  as  the  grid  is  repaired  after 
an  attack  (which  can  make  slow-to-repair  grid  components  more  attractive 
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targets).  An  attack  is  optimal  in  the  first  model  if  it  maximizes  the  total, 
instantaneous,  unserved  demand  for  power;  in  the  extended  model,  it  is  optimal  if 
it  maximizes  the  total  cost  of  unserved  demand  for  energy.  These  models  have 
been  tested  using  standard  reliability  test  systems  (RTSs)  drawn  from  IEEE 
(1999-1,  1999-11). 

This  thesis  makes  no  attempt  to  extend  the  modeling  approach  for  power 
flows  used  by  Salmeron  et  al.  Thus,  some  values  representing  power-system 
behavior  (such  as  reactive  power  flows,  line  losses,  voltage  magnitudes, 
transformer  tap  positions,  etc.)  are  still  disregarded  and/or  are  assumed  fixed  in 
this  work.  However,  we  assess  the  accuracy  of  that  approach  by  comparing 
results  from  our  approximating  “DC  power-flow  model”  with  those  provided  by  a 
full  AC  power-flow  model.  We  use  the  AC  model  embedded  in  the  PowerWorld 
Simulator  (PowerWorld  2003)  to  carry  out  this  comparison. 

Salmeron  et  al.  (2003)  propose  heuristic  algorithms  for  finding  an 
approximate  solution  to  their  interdiction  models.  These  heuristic  solutions  are 
generally  within  10%  of  optimality  but  in  some  instances  that  gap  can  reach  25%. 
The  possibility  of  large  gaps  was  discovered  after  publication  of  the  referenced 
paper,  when  the  authors  were  able  to  solve  equivalent  Mixed  Integer  Linear 
Program  (MIP)  formulations  exactly  (Salmeron  et  al.  2004).  These  MIP 
formulations  have  been  reviewed  and  consolidated  in  the  framework  of  this 
thesis 

In  theory,  a  MIP  can  be  solved  exactly  by  standard  techniques  such  as 
“branch-and-bound”  (e.g.,  Wolsey  1998,  pp.  95-107).  However,  these 
techniques  cannot  always  solve  the  large-scale  problems  found  in  practice,  and 
this  motivates  the  search  for  alternative  techniques  with  better  scalability. 
Benders  decomposition  (Benders  1962)  is  a  well-known  method  for  solving  large- 
scale  MIP  models,  and  this  thesis  investigates  its  use  to  solve  the  interdiction 
problem  efficiently.  Previous  application  of  Benders  decomposition  to  interdiction 
problems  can  be  found  in  Cormican  (1995)  and  Israeli  and  Wood  (2002). 
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B.  THESIS  OUTLINE 

Subsequent  chapters  in  this  thesis  are  organized  as  follows:  Chapter  II 
introduces  the  DC  power-flow  model  that  is  the  basis  for  the  interdiction  model.  It 
also  defines  the  interdiction  model  and  presents  the  heuristic  and  MIP 
approaches  developed  previous  to  this  work  to  solve  the  problem.  Finally, 
Chapter  II  introduces  Benders  decomposition  and  its  associated  algorithm  in  the 
framework  of  our  interdiction  problem.  Chapter  III  validates  the  DC  power-flow 
model  through  comparison  with  an  AC  power-flow  model.  We  compare  power 
flows  for  different  interdiction  plans  on  IEEE’s  “One  Area  RTS.”  Chapter  IV 
details  the  application  of  Benders  decomposition  to  our  interdiction  model.  A 
three-bus  example  demonstrates  the  decomposition  process,  which  is  then 
applied  to  the  interdiction  problem  without  system  restoration.  Computational 
results  are  presented.  Chapter  V  presents  the  Benders  decomposition  for  the 
interdiction  model  with  system  restoration  over  time.  Chapter  VI  discusses 
refinements  to  the  algorithm  and  provides  computational  results.  Chapter  VII 
summarizes  results  and  recommends  topics  for  future  work.  Appendices  A  and 
B  show  the  full  derivation  of  the  interdiction  models,  which  are  introduced  in 
Chapter  II  and  further  described  in  Chapters  IV  and  V.  Appendix  C  describes  the 
linearization  of  cross  products  for  the  three-bus  problem  developed  in  Chapter 
IV. 
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II.  PRELIMINARIES 


This  chapter  introduces  the  DC  Optimal  Power  Flow  model  (DCOPF)  that 
is  the  basis  for  the  development  of  two  interdiction  models:  without  system 
restoration  (see  Section  B),  and  including  system  restoration  (see  Section  C). 
Both  Sections  B  and  C  describe  heuristic  and  MIP  approaches  explored  prior  to 
this  research  for  solving  the  interdiction  models.  Section  D  introduces  Benders 
decomposition  in  the  context  of  our  models. 

A.  INTRODUCTION  TO  DCOPF 

At  the  heart  of  the  interdiction  models  is  a  DC  power-flow  model  (DCPF), 
so  we  first  provide  some  background  on  this  and  related  models.  DCPF 
simplifies  the  so-called  “full  AC  power-flow  model”  (ACPF),  which  is  generally 
accepted  as  a  valid  representation  of  power  flows  under  steady-state  conditions 
and  symmetry  (Frauendorfer  et  al.  1993).  The  DCPF  representation  entails 
various  assumptions  which  may  be  acceptable  in  the  context  of  security  analysis 
(Wood  and  Wollenberg  1996).  For  example,  DCPF  disregards  reactive  power 
effects  and  assumes  nominal  voltage  magnitudes  at  the  buses. 

DCPF  can  be  used  to  construct  an  optimization  problem,  DCOPF,  which 
optimizes  a  merit  function  subject  to  DCPF  constraints.  The  merit  function 
typically  measures  generating  costs.  In  contrast,  our  DCOPF  not  only  minimizes 
generation  costs,  but  also  the  penalty  associated  with  unmet  demand  (“load 
shed”),  because  we  cannot  guarantee  that  all  demand  will  be  met  in  the 
presence  of  interdiction.  The  DCOPF  formulation  from  Salmeron  et  al.  (2003)  is 
presented  below  for  completeness.  (Definitions  of  terms  can  be  found  in  the 
glossary.) 

Sets: 

ze/,  set  of  buses 

g  G  G ,  set  of  generating  units 
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/el,  set  of  transmission  lines 

c  e  C ,  set  of  consumer  sectors 

se  S ,  set  of  substations 

/e  ,  subset  of  buses  at  substation  5 

g&Gi,  subset  of  generating  units  connected  to  bus  i 

/e  if" ,  subset  of  lines  connected  to  bus  / 

subset  of  lines  connected  to  substation  5  (including 
transformers,  which  are  represented  by  lines) 

G*  c  G ,  <^L ,  I*  I ,  S*  <^S ,  interdictable  generators,  lines,  buses,  and 

substations,  respectively.  These  are  “interdictable 
components.” 

Parameters  (units,  if  applicable): 

o(/),  d{l),  origin  and  destination  buses  of  line  /  respectively  (more  than 
one  line  with  the  same  o(/) ,  d{l)  may  exist). 

i{g) ,  bus  for  generator g ,  i.e.,  ge 

s{i),  substation  s&S  associated  with  bus  /g  f 

d^^ ,  load  of  consumer  sector  c  at  bus  i  (MW) 

^Line  ^  transmission  capacity  for  line  I  (MW) 

,  maximum  output  from  generator  g  (MW) 

resistance,  and  reactance  of  line  /  respectively  (Q).  (We 
assume  x,  »  r,  ) 

Bi ,  series  susceptance,  calculated  as  B,  =  x,  lirf-  +  xf  ( 
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generation  cost  for  unit  g  ($/MWh) 


load-shedding  cost  for  customer  sector  c  at  bus  i  ($/MWh) 


Decision  variables  (units): 

pGen  ^  generation  from  unit  g  (MW) 

,  power  flow  on  line  /  (MW) 

,  load  shed  (unmet)  for  customer  sector  c  at  bus  i  (MW) 

,  phase  angle  at  bus  i  (radians) 


Remark:  All  units  are  converted  to  per  unit  (p.u.)  values  for  a  base  load  of 
100  MW.  (Remark:  We  use  this  same  value  in  all  implemented  models.)  For 
example,  to  convert  r,  to  per-unit  values,  we  consider  the  transmission  line  (/) 

nominal  voltage  rating  E,  (in  kV).  Then: 


r,  (per  unit)  =  (ohms)  x  100  (MW)/(line  rating(kV))^  = 


lOOr^ 

( MVA 

_  lOOr^ 

_  lOO'i  1 

Ef 

1  (kV)^  ; 

“  Ef 

“  Ef  ' 

1  V  J 

lOOr 

— ^p.u.  (since  lV=lAxlt2) 


The  same  process  is  carried  out  for  reactance  (x, )  conversions.  For  the 
rest  of  MW  magnitudes,  we  use  per-unit  values  after  dividing  by  the  100  MW 
base  load. 

The  formulation  of  DCOPF  is: 


(DCOPF):  min 

pGen  pLine  ^  q 


g  I  c 


s.t. 


V/e  L 

(2.1) 

^  pGen  ^  pLine 

X  Pr  = 

Vi 

(2.2) 

l\d(l)^i  c 

pLine  ^  pLine  ^  pLine 

yiGi 

(2.3) 
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\fi,c 


(2.4) 

(2.5) 


pGen  ^  pGen 
—g  ~  ^g 

0<5, 


<P. 


Gen 


This  model  assumes  fixed  (nominal)  voltage  magnitudes  and  does  not 
represent  shunts,  reactive  power  flows,  power  losses,  DC  lines,  transformer  tap 
positions,  phase  transformers,  and  other  features  that  could  destroy  the  linear 
structure  of  the  approximation.  Discrepancies  in  power  flows  in  DCOPF  and  a 
full  ACPF  are  discussed  in  Chapter  III. 


B.  THE  INTERDICTION  MODEL  WITHOUT  SYSTEM  RESTORATION 


1.  Interdiction  Model 

After  some  manipulation,  DCOPF  above  can  be  stated  in  standard  form 


(DCOPF):  min  cy 

y 

s.t.  Ay  =  b  (2-6) 

J>0 

where  y  represents  generation  outputs,  load  shedding,  power  flows  and  phase 
angles  variables,  along  with  slack  variables  associated  with  constraints  (2.3)- 
(2.5). 


The  interdiction  model  extends  (2.6)  in  order  to  choose  the  most 
disruptive,  resource-constrained  interdiction  plan,  ^gA,  where  A  is  a  set  of 
binary  vectors  representing  possible  terrorist  attack  plans.  An  interdiction  plan  is 
represented  by  the  binary  vector  S ,  where  is  1  if  component  e  of  the  power 

grid  is  attacked  and  is  0  otherwise.  We  assume  that  terrorists  attempt  to 
maximize  the  minimum  post-interdiction  “disruption”  which  is  measured  as  the 
generating  cost  plus  the  load-shedding  cost  (penalty)  evaluated  by  DCOPF. 


A  simple  representation  of  the  max-min  Interdiction  of  the  DCOPF  (I- 
DCOPF)  model  is  shown  below: 
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(I-DCOPF)  :  max  mine  y 

SeA  y 

s.t  giS,y)  =  0  (2.7) 

J>0 

where  g{S,y)  represents  a  set  of  constraints,  some  of  which  involve  functions 

that  are  nonlinear  in  iS,y)  (see  Appendix  A.1).  For  a  given  interdiction  plan,  S , 
the  inner  minimization  is  still  a  DCOPF  model  which  is,  of  course,  linear  in  y . 

Prior  to  this  thesis,  model  (2.7)  was  tackled  by  (a)  using  a  heuristic 
algorithm  and  by  (b)  reformulating  (2.7)  as  a  MIP  which  was  solved  with  a 
standard  solver.  This  thesis  contributes  to  (b)  by  refining  the  MIP  formulation 
that  serves  as  basis  for  the  Benders  decomposition  approach.  In  order  to  place 
this  thesis’  contributions  in  context,  the  following  two  subsections  summarize 
approaches  (a)  and  (b). 

2.  Heuristic  Approach 

A  heuristic  algorithm  introduced  by  Salmeron  et  al.  (2003)  solves  model 
(2.7)  approximately,  and  has  been  tested  using  small-  and  medium-scale 
networks  drawn  from  the  IEEE  Reliability  Test  Data  (1999-1,  1999-11).  This 
heuristic  solves  DCOPF  for  a  given  grid  configuration  (i.e.,  for  a  given  interdiction 
plan).  Then  it  assigns  “values”  to  interdictable  electrical  components  in  the  grid 
(lines,  buses,  substations,  etc.)  based  on  present  and  previous  flow  patterns 
associated  with  the  components.  Finally,  the  heuristic  maximizes  the  value  of 
the  components  to  be  interdicted  while  excluding  previously  explored  solutions. 
This  process  is  repeated  for  a  pre-specified  number  of  iterations.  This 
decomposition-based  approach  has  been  empirically  proven  to  obtain  good,  but 
not  necessarily  optimal,  interdiction  plans.  Figure  1  depicts  the  basic  cycle  of  the 
algorithm.  Additional  details  on  the  heuristic  can  be  found  in  the  cited  reference. 
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Figure  1 .  Framework  for  a  heuristic  interdiction  algorithm.  (From  Salmeron  et  al. 

2003.) 

The  inability  to  prove  optimality  of  the  heuristic  solution  (unless  all  possible 
solutions  are  enumerated  explicitly)  has  led  to  the  search  for  alternative 
representations  of  (2.7)  that  are  amenable  to  solution  by  MIP  techniques,  such 
as  branch-and-bound,  as  described  in  the  following  section. 

3.  The  MIP  Approach 

a.  Overview 

The  development  of  a  MIP  model  enables  validation  of  the  heuristic 
solution  and  the  development  of  exact  decomposition  algorithms  that  are 
frequently  used  solve  large-scale  problems.  (Typically  the  solutions  are 
approximate,  but  then  large-scale  systems  are  almost  never  solved  exactly  by 
any  technique.)  The  derivation  of  the  MIP  begins  with  the  linearization  of 
constraints  in  g{S,y)  =  0,  followed  by  the  dualization  of  the  inner  minimization 
problem  in  (2.7).  This  procedure  converts  the  max-min  problem  into  a  nonlinear 
maximization  problem  which  is  linearized  through  additional  steps.  Both  of  these 
linearization  techniques  are  fully  explained  in  Salmeron  et  al.  (2004),  but  are 
briefly  described  in  the  following  for  continuity. 

b.  Derivation 

The  nonlinearities  in  the  expression  of  g{S,y)  =  0  in  model  (2.7)  are 
associated  with  admittance  equations  in  the  presence  of  interdiction  variables 
These  have  the  form: 


10 


no -•5") 

'y  site  4  J 


(2.8) 


Note  that  in  the  interdiction  model,  the  original  admittance  equation 
(DC.1)  for  every  interdictable  line  /  is  modified  as  (2.8)  to  account  for  potential 
interdiction  (directly  or  indirectly)  of  the  line,  in  which  case  the  admittance 
equation  should  not  be  enforced.  This  leads  to  nonlinearities  in  the  form  of 
multiple  cross-products. 


Following  Salmeron  et  al.  (2004),  (2.8)  can  be  fully  linearized.  For 
example,  one  can  convert  the  nonlinear  constraint  P  =  B{0^  - - d^){\  - 5^)  into 
the  following  two  linear  inequalities: 

P-B{O^-0,)<M{S,+S,) 

P-B(0^-0,)>-M(S,+S,) 

Here,  M  =  P  +  B0  ,  where  0  is  an  upper  bound  on  the  maximum 
phase  angle  difference  between  buses  a  and  b.  (2.9)  enforces  P  =  B{0^-0J 
when  all  the  related  (^-variables  are  0,  and  eliminates  the  constraint  when  any  of 
these  S  variables  is  1.  The  resulting  model  is  still  called  l-DCOPF  in  this  thesis, 
and  is  specified  in  detail  in  Appendix  A.1 . 


We  proceed  by  taking  the  dual  of  the  inner  minimization  in  (2.7), 
assuming  constraints  in  g{S,y)  =  0  have  been  linearized  as  described  above. 
This  converts  the  max-min  problem  into  a  simple  maximization.  The  resulting 
model  is  called  DI-DCOPF  (“Dual  of  l-DCOPF,”  although  only  the  inner 
minimization  is  dualized).  The  full  formulation  is  given  in  Appendix  A.2.  This 
model’s  objective  function  consists  of  a  linear  function  of  continuous  vector  tt, 
along  with  nonlinear  costs  involving  cross-products  of  type  Stt.  The  new 
variables  correspond  to  duals  of  each  of  the  balance  and  bound  constraints  in 
l-DCOPF.  The  terms  associated  with  the  Stt  products  arise  from  the  line  and 
generator  capacity  constraints. 
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The  resulting  problem  can  be  stated,  in  brief,  as: 

(DI-DCOPF):  maxmax  b^7r+b''v 

K 

s.t.  K7C<c  (y)  (2.10) 

VG  V{S,:7r) 

Here,  Zj'^and  F  are  appropriately  dimensioned  coefficient  vectors, 
;r  is  a  vector  of  dual  variables  to  constraints  in  l-DCOPF  and  v  is  a  vector  where 
each  component  is  the  result  of  a  cross  product  of  a  variable  S  and  a  variable 
TT .  Variables  in  parentheses,  here  and  elsewhere,  are  dual  variables  for  the 
associated  constraints;  in  this  case  the  variables  are  y . 

To  linearize  a  cross-product  involving  a  variable  and  a 
;r-variable,  we  create  an  intermediate  variable,  e.g.,  =  with  the 

property  that,  Sj  =  0^Vj^=0,  and  =  1  ^  =;r^ .  Since  the  sign  of  is 
known  (assume  ^^>0  for  demonstration  purposes),  we  add  the  following  linear 
constraints,  represented  as  vGV(S,7r),  to  ensure  the  above  relationships 
between  and  hold: 

Vjk  ^  ^k 

Vjk  =  ^J^k  =  '  Vjk  ^  ^k-^k(^-^j)  (2-1 1 ) 

0<7rk<7rk 

Vjk^O, 

where  is  an  appropriate  upper  bound  on  . 

Let  us  drop  the  {j,k  )  subindices  for  simplicity.  Then,  the  following 
scheme  depicts  the  validity  of  the  representation  in  (2.1 1 ): 
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v<M 

v<n 

v>7t-7i{V-5) 

^<7C<7i 

v>0 


Remark:  When  ;r<  0 ,  the  linearization  is  carried  out  as: 

Vjk  ^ 

Vjk=^j^k  =  ^^jk^^k+^k(^-^j)  (2-13) 

-^,<7r,<0 

Vjk  ^  0 

Finally,  generation  phase  angles,  power  flows  and  load  shedding 
are  jointly  denoted  by  the  dual  variable  y . 

The  dualization  and  linearization  described  above  allow  us  to 
represent  the  interdiction  problem  as  a  MIP.  The  resulting  model  is  called  LDI- 
DCOPF  (Linearization  of  DI-DCOPF;  see  Appendix  A. 3). 

C.  INTERDICTION  MODEL  WITH  SYSTEM  RESTORATION 

1.  System  Restoration 

The  model  described  in  Section  B  is  a  static  representation  of  power 
disruptions  at  a  given  point  in  time.  That  model  does  not  represent  the 
consequences  of  the  variability  in  repair  times  of  damaged  system  components 
and  the  increased  cost  of  load  shedding  resulting  from  repair-time  delays. 
Unless  the  outage  duration  of  each  interdictable  component  is  the  same,  the 
interdiction  model  presented  in  Section  B  cannot  capture  this  effect.  For 
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example,  Salmeron  et  al.  (2003)  show  that  if  the  interdictor’s  decisions  are  driven 
by  the  goal  of  maximizing  short-term  power  shedding  (or  its  hourly  cost),  the 
resulting  interdiction  plan  might  be  far  from  optimal  had  his  or  her  goal  been 
maximizing  total  energy  shed  (or  its  cost)  until  the  system  repair  is  completed. 
This  issue  is  especially  important  when  transformers  and  other  difficult-to-replace 
equipment  are  subject  to  interdiction. 

Figure  2  depicts  of  the  ideas  above  using  two  graphs. 


Figure  2.  Power  disruption  and  energy  disruption.  The  model  without  restoration 
provides  the  optimal  interdiction  plan  according  to  instantaneous  picture  of 
power  shed  (left).  The  model  with  restoration  over  time  accounts  for 
energy  disruption  (right).  (From  Salmeron  et  al.  2004.) 


2.  Interdiction  Model 

Salmeron  et  al.  (2004)  show  how  l-DCOPF  (2.7)  can  be  extended  to 
handle  the  time  associated  with  repairs.  This  resulting  model  seeks  to  maximize 
total  cost  of  energy  shed,  accounting  for  the  fact  that  the  system  will  be 
progressively  restored  over  time.  The  interdiction  model  with  system  restoration 
can  be  stated  in  brief  as: 
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max  min  }  D,  cy, 

Se^  y  ^ 

sX.  gXS,y,)  =  0,  yt&T  (2.14) 

J,>0,  V?er 

This  model  uses  interdiction  constructs  to  couple  instances  of  DCOPF, 
one  for  each  system  state  that  represents  a  stage  or  “time-period”  of  system 
repair:  Note  the  time  subindex  t,  where  T  is  the  set  of  time  periods  based  on 
repair  times  of  all  components  in  the  system  (see  Appendix  B.1).  That  index  is 
added  to  all  of  the  inner  decision  variables  to  account  for  power  flows  in  the  grid 
and  other  decision  variables  that  change  over  time.  represents  the  duration 

of  time  period  t,  y,  is  the  same  y  described  in  (2.6)  but  for  every  time  period  t, 
and  g,  is  the  resulting  set  of  nonlinear  equations  in  {S,y,).  Note  the  time- 
dependent  constraints  g,{S,y,)  =  0  not  only  ensure  that  some  components  are 
out  of  service  right  after  interdiction,  but  also  guarantee  that  components  will  be 
in  service  again  after  their  posted  repair  time.  Thus,  the  constraints  g,(^,j,)  =  0 
for  different  periods  t  are  coupled  by  the  variables  S .  The  full  model  derivation, 
described  in  detail  in  Salmeron  et  al.  (2004),  consists  of  the  inner  DCOPF  model 
with  system  restoration  (referred  to  as  DCOPF-R)  and  the  associated  interdiction 
model  (referred  to  as  l-DCOPF-R);  see  Appendix  B.1  for  additional  definitions 
and  notation,  and  Appendix  B.2  for  the  formulation. 

As  in  the  case  without  system  restoration,  a  solution  to  the  interdiction 
model  with  system  restoration  can  be  obtained  through  a  heuristic  algorithm  or 
through  a  MIP  reformulation  of  the  problem,  as  described  in  the  following 
paragraphs. 

3.  Heuristic  Approach 
a.  Description 

The  Heuristic  for  a  model  with  system  restoration  over  time 
(schematically  depicted  in  Figure  3)  follows  the  approach  of  that  for  a  model 
without  system  restoration. 
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For  a  fixed  interdiction  vector  S  (i.e.,  S  =  S),  the  inner  DCOPF-R 
subproblem  provides  the  joint  flow  patterns  for  a  number  of  system  stages 
(restoration  periods).  The  subproblem  decomposes  into  |r|  sub-subproblems 
(because  when  S  is  fixed,  the  subproblem  is  no  longer  coupled  by  time  periods). 
Each  subproblem  consists  of  an  instance  of  DCOPF  with  some  subset  of  system 
components  being  “out  of  service.”  Components  that  are  out  of  service  are 
determined  “on  the  fly”  as  the  algorithm  solves  a  DCOPF  model  for  every  time 
period. 

Salmeron  et  al.  (2004)  redefine  the  concept  of  a  component’s  value 
in  this  dynamic  model  by  factoring  the  time  it  would  be  out  of  service  if 
interdicted. 


r 


Solve  the  DC-OPF-R  for  the  given 
Interdiction  Plan: 

Solve  \T\  DC-OPF  problems 


Based  on  present  and  previous  flow 
patterns,  assign  an  (energy-based) 
“Value”  to  each  interdictable  asset 


MP-R:  Maximize  the  Value  of  the 
Assets  to  be  Interdicted  (excluding 
previously  explored  solutions) 


Figure  3.  Interdiction  algorithm  framework  with  restoration.  DCOPF-R  decomposes 
into  |r|  sub-subproblems,  each  being  an  instance  of  DCOPF  with  some 
components  “out  of  service.”  (From  Salmeron  et  al.  (2004).) 


4.  The  MIP  Model 
a.  Overview 

The  structure  of  l-DCOPF-R  allows  us  to  convert  it  into  a  standard 
MIP  through  a  process  of  dualization  and  linearization,  just  as  in  the  case  without 
system  restoration. 
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b.  Derivation 

As  the  new  interdiction  model  is  a  variant  of  the  model  without 
system  restoration,  where  the  inner  model  has  a  similar  structure  to  a  power  flow 
model  (replicated  over  time),  the  development  of  the  MIP  formulation  is 
performed  as  in  Section  B:  First,  admittance  equations  are  linearized  (Appendix 
B.2).  Next,  the  inner  minimization  problem  is  dualized  (Appendix  B.3).  Finally, 
cross-products  in  the  objective  function  are  linearized.  The  resulting  MIP,  called 
LDI-DCOPF-R,  is  specified  in  Appendix  B.4. 

D.  BENDERS  DECOMPOSITION 

Benders  Decomposition  is  a  well-known  technique  for  solving  MIPs 
(Benders  1962).  This  technique  has  proven  effective  in  solving  large-scale 
optimization  problems  (e.g.,  facility  location  problems,  two-stage  stochastic 
optimization  problems,  etc.)  including  interdiction  problems  (Cormican  1995, 
Israeli  and  Wood  2002).  Benders  decomposition  is  well  suited  when  the 
assignment  of  a  subset  of  variables,  so-called  “complicating  variables,”  yields 
easily  solvable,  disconnected  and  convex  subproblems,  such  as  the  network-flow 
problems  in  Geoffrion  and  Graves  (1974)  and  Brown  and  McBride  (1984);  these 
MIPs  could  probably  not  have  been  solved  directly,  i.e.,  using  branch  and  bound. 
Based  on  the  successful  application  of  Benders  decomposition  to  these  similar 
problems,  we  implement  it  to  solve  LDI-DCOPF  and  LDI-DCOPF-R. 

Since  the  development  of  Benders  decomposition  is  similar  for  both 
models,  we  focus  on  LDI-DCOPF.  The  MIP  formulation  of  this  model  embeds  a 
bi-level  structure  that  makes  Benders  decomposition  suitable  for  solving  the 
problem;  The  first  level  is  the  interdictors’  level,  which  is  modeled  through  binary 
decision  variables.  The  second  level  is  an  instance  of  DCOPF  (actually,  its  dual), 
where  some  components  of  the  original  power  grid  have  been  removed  through 
the  first-level  decision  variables.  Because  DCOPF  is  linear,  it  satisfies  Benders’ 
theoretical  requirement  for  convexity  in  the  subproblem. 
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Our  interdiction  problem  (LDI-DCOPF)  may  be  stated,  in  brief,  as: 

(P):  z*  =  max  c'’v  +  c^Tt 

S,v,7r 

s.t.  A^v^A^tc^  B5  =  b  {y) 

;r>0  (2.15) 

v>0 

Sg  A 

Note:  This  problem  is  referred  to  as  LDI-DCOPF  elsewhere  in  this  thesis, 
but  for  simplicity  in  this  section,  we  refer  to  it  as  problem  (P). 

Note:  non-positive  tt's  and  v's  have  been  converted  to  non-negative  tt's  and 
v's . 

In  our  Benders  decomposition  scheme,  the  interdiction  vector  S  plays  the 
role  of  the  vector  of  “complicating  variables.” 

Rewriting  (P)  as: 

(P):  max  max  c‘’v-i-c'^;r 

SeA  v,7r 

s.t.  A'’v  +  A^n:  =  b  -  BS  {y)  ^2  16) 

v>0 
;r>  0 

the  dual  of  the  inner  maximization  above  is: 

min  y{b-BS) 

y 

s.t.  yA^>c^  (v)  (2.17) 

(TT) 

where  y  represents  the  vector  of  dual  variables  corresponding  to  the  constraints 
of  the  primal  problem  in  (2.16).  We  assume  the  polytope 
H  =  [y\yA'' ^cyyA’"  always  contains  a  feasible  solution,  and  its  optimal 
solution  is  not  unbounded.  We  will  show  later  that  these  assumptions  are  valid. 

Since  the  solution  to  (2.17)  will  lie  at  one  of  the  extreme  points  of  the 
feasible  region  H ,  we  can  rewrite  (2.17)  as: 
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mmy^{b- B5)  (2.18) 

iel 

where  I  =  {\,2,...,k}  indexes  the  set  of  extreme  points  of  H  which  is 


Since  mmy-{b-Bd)  <  y.(b-BS),  I ,  we  can  rewrite  (P)  as: 

iel 


(P) :  max  z 

6gA,zgE 


(2.19) 


s.t.  z  <  yi{b - Bd)  V/ e  / 

If  we  relax  (2.19)  by  considering  a  subset  of  constraints  zg/,  indexing 
these  constraints  by  z  =  1,2,...,A: ,  at  a  given  iteration  k  of  our  Benders 
Decomposition  Algorithm  (BDA),  we  have  the  following  “master  problem”  (MP)  to 
solve: 


(MP,):  z,  =  maxz 

6gA,zgIR 

s.t.  z  <  y-{b-BS) 


i  =  h2, 


k 


(2.20) 


(MP^)  is  a  relaxation  of  (P),  and  therefore  its  optimal  objective  function 


value,  z^ ,  constitutes  an  upper  bound  on  z* ,  the  optimal  objective  value  to  (P). 


If  the  solution  to  (MP^),  denoted  (4’4)>  'S  not  optimal  to  (P),  then 
must  violate  some  constraints  in  (2.19)  not  included  in  (MPJ,  i.e.,  constraints 
associated  with  extreme  points  }  ■  (Remark:  Here,  the  subindex  k 

refers  to  the  iteration  counter,  and  is  the  incumbent  interdiction  vector.) 

A  natural  way  to  check  for  one  of  these  extreme  points  is  to  identify  the 
vector  satisfying  that  violates  z^  <  the  most.  That 

identification  problem  can  be  stated  as  the  following  optimization  “subproblem:” 

(SP,(4)):  mmy,^,{b-Bd,) 

rt+i 

sX.  y,,,A^>c^  (V,)  (2.21) 

{7U,) 


or,  in  its  primal  version: 
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(2.22) 


(SP^(4,)):  max  c'’v  +  c’^7U 

V,7I 

s.t.  A'v-^  =  b  -  bSj^  (Pi:+i) 

v>0 
;r>  0 

It  is  interesting  to  note  that  model  (SPJ  in  (2.22)  is  the  same  problem  as 
problem  (P)  in  (2.16)  when  the  value  of  the  decision  variable  S  is  fixed  as  S  =  d,^. 

(SPJ  is  therefore  feasible,  because  problem  (P)  is  feasible  for  any  S  =  Si^. 

(Recall  that  the  inner  problem  in  (2.16)  is  equivalent  to  a  DCOPF  model,  which  is 
always  feasible,  even  if  this  feasibility  entails  penalties  for  unmet  load.)  In 

addition,  (SP^)  is  bounded  for  any  ^  =  4-  because  our  DCOPF  models  have 
non-negative  cost  coefficients  only. 

The  solution  to  (SPJ,  denoted  (v^,i^t),  along  with  4  from  model  (MP^) , 
form  a  combined  solution  (4,v^,#J,  which  is  feasible  to  our  original  (P).  Its 
objective  function  value  is  c'’v^+c''Jr,^,  which  represents  a  lower  bound  on  the 
optimal  solution  to  problem  (P). 

Solving  the  subproblem  above  yields  a  new  extreme  point  whose 
associated  cut  z<  will  be  violated  by  the  incumbent  solution  (4^4) 

unless  (4^4)  is  already  optimal  to  (P).  The  newly  generated  cut  is  added  to 
(MPJ ,  becoming  (MP^^j),  and  the  process  is  repeated. 

E.  DECOMPOSITION  ALGORITHM 

The  mathematical  derivation  described  in  Section  D,  leads  to  the  following 
step-by-step  algorithm:  (The  r?()  notation  refers  to  the  optimal  objective 
function  value  of  the  problem  in  the  argument.) 

Benders  Decomposition  Algorithm  (BDA): 

Input:  Initial  solution  4^^>  matrices  A\  A’' ,  and  B,  vector  b,  and  an 
optimality  tolerance  f>0. 
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Output:  f -optimal  interdiction  plan  S* ,  and  associated  power  flows, 
generation  phase  angles,  and  load  shedding  (jointly  denoted  as  /). 

1.  Set  the  iteration  counter  k\  =0,  the  lower  bound  LB:=-oo,  and  the 
upper  bound  UB:=-i-oo. 

2.  Solve  SP.(4)  for  (V., Let  LB.  :=i>(sp,(4)) . 

3.  If  LB^<LB,  then  update  the  lower  bound:  LB;=LB^,  and  set 
5  :=d„y  :=y,^,. 

4.  If  UB  -  LB  <  f ,  STOP ;  otherwise  continue  to  step  5. 

5.  Add  the  newly  generated  cut:  z  <  y^^^(b - to  (MP^). 

6.  Set  k:=k  +  \. 

7.  Solve  (MPJ  for  4  >  and  4  =  ^(MPJ  . 

8.  Update  the  upper  bound:  UB  :=  ^ ,  and  return  to  step  2. 

At  each  iteration,  the  solution  to  SP(4)  gives  us  a  new  candidate 

solution  (4’A+i)-  Whenever  the  objective  value  for  this  candidate  solution  is 
greater  than  the  previous  lower  bound,  we  replace  the  lower  bound  with  this 
value.  This  provides  a  non-decreasing  lower  bound.  The  solution  to  (MP^)  is 

monotonically  decreasing  with  each  cut  added,  and  therefore  we  automatically 
update  the  upper  bound  at  each  iteration.  This  procedures  is  repeated  until  the 
bounds  converge,  or  are  sufficiently  close,  which  in  the  worst  case  will  occur  after 
generating  all  possible  extreme  points  of  the  subproblem,  i.e.,  when  k  =  K . 

To  improve  the  BDA  efficiency,  we  incorporate  some  enhancements  in 
Chapter  VI.  Figure  4  depicts  the  original  BDA  for  ease  of  comparison  with  the 
flow  chart  for  the  enhanced  BDA. 
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Benders  Decomposition  Aigorithm  (BDA) 


LB:=-oo 

UB:=oo 

k:=0 


Figure  4. 


Benders  Decomposition  Algorithm  flowchart, 
22 


Using  the  above  template,  we  implement  BDA  using  the  full  formulation  of 
the  interdiction  problem,  with  and  without  system  restoration,  and  test  it  on 
several  benchmark  cases. 

The  decomposition  for  the  model  without  system  restoration  is  presented 
in  Chapter  IV,  starting  with  a  simple  example  to  familiarize  the  reader  with  the 
process.  Chapter  V  expands  upon  this  model  to  take  into  account  system 
restoration,  and  presents  computational  results  for  that  model  for  BDA. 

Before  describing  details  of  the  decomposition  process,  we  validate  the 
inner  power-flow  model  in  the  next  chapter. 


23 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


24 


III.  VALIDATING  THE  DC  POWER  FLOW  MODEL 


This  chapter  validates  the  DC  Power  Flow  model  (DCPF)  introduced  in 
Chapter  II,  whose  accuracy  must  be  established  to  ensure  the  reasonableness  of 
our  interdiction  model,  l-DCOPF.  We  do  this  using  non-interdicted  and 
interdicted  instances  of  the  IEEE  One  Area  test  case  (IEEE  Reliability  Test  Data 
1999-1). 

A.  INTRODUCTION 

DCPF  is  a  simplification  of  the  full  AC  Power  Flow  model  (ACPF)  (e.g.. 
Wood  and  Wollenberg  1996).  ACPF  includes  reactive  power  flows  (DCPF 
disregards  these),  voltage  magnitudes  at  the  buses  (DCPF  assumes  nominal 
voltages,  i.e.,  1.00  p.u.),  and  power  losses  on  the  line  (DCPF  assumes  these  to 
be  zero),  among  others. 

Overbye  et  al.  (2004)  conduct  an  analysis  to  validate  the  use  of  DCPF  in 
place  of  ACPF  for  market  analyses;  they  find  it  adequate.  Wood  and  Wollenberg 
(1996)  indicate  that  DCPF  is  adequate  for  security  analysis  of  many  systems. 
Although  we  believe  that  similar  results  should  hold  for  interdiction  problems, 
further  analysis  is  desired  because  these  problems  have  additional 
complications.  In  particular,  they  involve  removing  critical  components  of  grids 
and  adjusting  loads  to  minimize  disruption  costs. 

We  validate  DCPF  through  comparative  analysis  with  the  ACPF  provided 
by  software  in  the  PowerWorld  Simulator  (PWS)  (PowerWorld  2003).  PWS  is 
standard  software,  widely  used  in  the  electric  power  industry.  Our  assessment  is 
carried  out  by  comparing  “optimal  power  flows”  (OPFs)  obtained  using  DCOPF 
with  those  identified  by  the  ACPF.  It  is  important  to  note  that  ACPF  does  not,  per 
se,  attempt  to  minimize  load  shedding  (or  its  cost)  as  DCOPF  does.  However, 
for  given  loads  and  generation,  it  provides  accurate  active  and  reactive  power 
flows  on  the  lines  and  bus  voltages  (magnitude  and  phase),  along  with  other 
system  values.  We  want  to  assess  whether  power  flows  and  voltage  angles. 
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obtained  through  DCOPF,  are  in  fact  close  to  those  from  ACPF  when  we  use  our 
optimized  generation  and  load  shedding  as  input  data. 

We  present  the  comparison  between  both  approaches  assuming  no 
interdiction  first,  followed  by  other  cases  where  critical  components  have  been 
interdicted. 


B.  DCOPF  VS.  ACPF  IN  THE  ABSENCE  OF  INTERDICTION 

To  enable  comparison  between  the  DCOPF  and  the  PWS  ACPF,  we 
incorporate  the  IEEE  RTS  One  Area  case  into  PWS  through  its  one-line  diagram 
building  interface.  First,  we  reconstruct  the  underlying  power  grid,  as  shown  in 
Figure  5,  using  the  following  symbols  to  represent  the  different  system 
components: 


Bus 


Line 


Fraction  of  line 

Fraction  of  line 

Generator 

capacity  in  use 

capacity  in  use 

(if  <80%) 

(if  >80%) 

Load 


vN/\/ 

/^vlv\ 


Transformer 


Flow  direction 


Circuit  breaker 
(close) 


Circuit  breaker 
(open) 


For  simplicity  (and  given  that  generating  units  are  not  interdictable  in  our 
test  cases),  we  aggregate  all  generating  units  at  a  bus  as  a  single  generator  at 
that  bus.  Next,  we  associate  grid  data  with  the  power  grid.  The  reconstructed 
one-line  diagram  includes  data  for  bus  nominal  voltages,  line  resistances, 
reactances  and  thermal  limits  (line  capacities),  aggregated  generating  capacities, 
etc.  Then,  we  incorporate  active  power  output  for  each  generating  unit  and 
active  load  met  at  each  bus  from  DCOPF  as  system  data  for  PWS’s  ACPF. 
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Finally,  we  run  ACPF  to  compute  actual  power  flows,  phase  angles  and  adjusted 
generation  at  the  swing  bus.  The  resulting  power  flows  on  the  lines  are 
compared  to  those  provided  by  DCOPF. 


The  one-line  diagram  in  Figure  5  represents  a  case  without  interdiction, 
which  we  call  “scenario  0,”  because  it  is  equivalent  to  a  case  with  zero 
interdiction  resource.  Section  C  analyzes  cases  with  interdiction. 


Figure  5.  IEEE  RTS  One  Area  case  one-line  diagram  without  interdiction.  This 

representation  is  constructed  using  the  PWS  software  (PowerWorld  2003). 
It  depicts  flow  directions,  fractions  of  line  capacities  used,  generation 

outputs  and  loads. 

A  side-by-side  comparison  of  the  resulting  flows  is  shown  in  Table  1.  We 
analyze  absolute  deviations  for  all  lines  and  transformers.  The  column  headed 
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“Line”  displays  the  line  name.  The  “From  Bus”  and  “To  Bus”  columns  indicate 
the  two  buses  connected  by  the  line.  The  “Transformer”  column  indicates  if  we 
are  using  a  line  to  represent  a  transformer.  The  “PWS  flow”  and  “PWS  loss” 
columns  show  active  power  flows  and  losses,  respectively,  obtained  with  PWS, 
while  the  “DCOPF  flow”  column  displays  active  power  flows  obtained  with 
DCOPF.  The  last  column  shows  the  percentage  absolute  deviation  between  the 
PWS  and  DCOPF  flows.  The  results  show  a  maximum  deviation  of  5.29%. 
Mean  deviation  is  1.18%  with  a  percentage  standard  deviation  of  0.012  %. 

The  example  above  indicates  that  DCOPF  yields  a  good  approximation  for 
the  true  AC  power  flow  in  this  RTS  case,  in  the  absence  of  interdiction.  In  the 
following  section,  we  analyze  what  happens  when  critical  elements  of  the  grid  are 
removed  through  interdiction. 


Line 

From 

Bus 

To  Bus 

Transformer 

PWS 

flow 

(MW) 

PWS 

loss 

(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A1 

101 

102 

No 

16.70 

0.01 

17.00 

1.80% 

A2 

101 

103 

No 

-26.50 

0.42 

-26.73 

0.87% 

A3 

101 

105 

No 

54.00 

0.68 

53.73 

0.50% 

102 

104 

No 

28.60 

0.29 

28.42 

0.63% 

102 

106 

No 

44.00 

1.01 

43.58 

0.95% 

103 

109 

No 

40.10 

0.53 

39.85 

0.62% 

103 

124 

Yes 

-244.30 

1.21 

-246.58 

0.93% 

104 

109 

No 

-45.30 

0.60 

-45.58 

0.62% 

105 

110 

No 

-17.20 

0.07 

-17.27 

0.41% 

106 

110 

No 

-96.80 

0.28 

-92.42 

4.52% 

A11 

107 

108 

No 

-123.50 

2.67 

-125.00 

1.21% 

108 

109 

No 

-151.90 

11.66 

-159.94 

5.29% 

108 

110 

No 

-130.50 

8.45 

-136.06 

4.26% 

A14 

109 

111 

Yes 

-153.80 

0.48 

-154.38 

0.38% 

109 

112 

Yes 

-185.20 

0.69 

-186.28 

0.58% 

110 

111 

Yes 

-203.10 

0.84 

-204.43 

0.65% 

110 

112 

Yes 

-234.20 

1.11 

-236.33 

0.91% 

111 

113 

No 

-223.90 

3.11 

-226.01 

0.94% 

111 

114 

No 

-132.20 

0.89 

-132.81 

0.46% 

112 

113 

No 

-170.00 

1.78 

-171.01 

0.59% 

112 

123 

No 

-245.20 

7.68 

-251.60 

2.61% 

113 

123 

No 

-183.20 

3.86 

-186.01 

1.53% 

114 

116 

No 

-322.20 

5.36 

-326.81 

1 .43% 

115 

116 

No 

45.60 

0.04 

45.90 

0.66% 

115 

121 

No 

-225.00 

3.14 

-227.24 

1.00% 
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Line 

From 

Bus 

To  Bus 

Transformer 

PWS 

flow 

(MW) 

PWS 

loss 

(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A26 

115 

124 

No 

248.10 

4.33 

246.50 

0.64% 

A27 

116 

117 

No 

-310.20 

2.96 

-312.52 

0.75% 

A28 

116 

119 

No 

87.40 

0.23 

86.81 

0.68% 

A29 

117 

118 

No 

-173.20 

0.61 

-172.69 

0.29% 

A30 

117 

122 

No 

-137.90 

2.78 

-139.82 

1.39% 

A31-1 

118 

121 

No 

-52.30 

0.08 

-52.85 

1.05% 

A32-1 

119 

120 

No 

-47.20 

0.11 

-47.19 

0.02% 

A33-1 

120 

123 

No 

-111.20 

0.38 

-111.19 

0.01% 

A34 

121 

122 

No 

-158.70 

2.35 

-160.18 

0.93% 

Table  1 .  ACPF  versus  DCOPF;  IEEE  RTS  One  Area  case  (no  interdiction).  For 
each  line  and  transformer,  we  show  the  PWS’s  ACPF  power  flow  and 
losses.  We  compare  power  flow  absolute  deviations  from  those  of  our 
DCOPF  model.  The  maximum  percentage  deviation  is  5.29%  for  the 
transmission  line  connecting  buses  108  and  109.  The  average  deviation 

is  1.18%. 


C.  DCOPF  VS.  ACPF  AFTER  INTERDICTION 

In  order  to  further  validate  DCOPF,  we  also  compare  flows  after 
interdiction.  The  main  effect  of  interdiction  is,  normally,  load  shedding.  We 
conduct  the  analysis  for  a  number  of  possible  scenarios;  recall  that  “scenario” 
corresponds  to  the  amount  of  interdiction  resource  available.  The  flow- 
comparison  procedure  is  repeated  for  the  RTS  One  Area  case  for  scenarios 
M  =  2,  4,  and  6 ,  in  l-DCOPF  shown  in  Appendix  A.1 . 

Although  overall  results  are  presented  for  all  cases,  we  will  illustrate  the 
one-line  diagram  and  comparison  table  for  scenario  M  =  6:  We  open  the 
interdicted  lines  (found  for  the  optimal  interdiction  plan  which  corresponds  to  this 
scenario)  by  opening  the  circuit  breakers,  at  the  ends  of  the  lines,  which  are 
represented  by  green  open  squares  in  Figure  6.  Then,  we  run  PWS  to  realize 
that  the  flow  in  a  number  of  lines  greatly  exceeds  the  lines’  capacities  (red  circles 
in  Figure  6).  While  a  line  can  exceed  its  nominal  capacity  temporarily,  the  long¬ 
term  emergency  rating  on  a  line  is  typically  20%  of  its  nominal  rated  capacity, 
and  a  line  can  handle  this  excess  flow  only  for  24  hours  (IEEE  1999-1). 
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U3  MW  -SSI 


Figure  6.  Effects  of  interdiction  to  IEEE  RTS  One  Area  case.  Some  lines  become 
overloaded,  which  is  unacceptable  but  for  a  few  hours  and  up  to  certain 
limits.  Note;  All  loads  are  kept  at  their  original  values  to  illustrate  the 
infeasibility  of  the  problem  after  interdiction. 

Having  determined  that  ACPF  shows  the  situation  to  be  infeasible  after 
interdiction,  we  reduce  the  loads  in  ACPF  to  match  those  provided  by  DCOPF 
and  run  PWS.  We  then  carry  out  a  comparison  of  the  flows  across  the  non- 
interdicted  elements  of  the  grid  with  those  provided  by  DCOPF;  see  Table  2. 

Note  that,  again,  all  the  DCOPF  flows  are  relatively  close  to  those 
produced  by  ACPF.  The  average  absolute  deviation  is  3.67%  with  a  standard 
deviation  of  0.0075%. 
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Line 

From  Bus 

To  Bus 

Status 

Transformer 

PWS 

flow(MW) 

PWS 

loss(MW) 

DCOPF 

flow(MW) 

Abs. 

deviation 

A1 

101 

102 

Closed 

No 

-3.60 

0.41 

-3.23 

10.28% 

A2 

101 

103 

Closed 

No 

68.30 

0.02 

66.01 

3.35% 

A3 

101 

105 

Closed 

No 

136.20 

1.18 

129.22 

5.12% 

A4 

102 

104 

Closed 

No 

104.80 

1.93 

99.95 

4.63% 

A5 

102 

106 

Closed 

No 

92.60 

9.82 

88.82 

4.08% 

A6 

103 

109 

Closed 

No 

76.70 

0.00 

76.52 

0.23% 

A7 

103 

124 

Closed 

Yes 

-10.40 

0.00 

-10.51 

1.05% 

A8 

104 

109 

Closed 

No 

95.70 

0.00 

99.95 

4.25% 

A9 

105 

110 

Closed 

No 

124.90 

0.00 

129.22 

3.34% 

A10 

106 

110 

Closed 

No 

82.80 

0.00 

88.82 

6.78% 

A11 

107 

108 

Open 

No 

0.00 

5.00 

0 

0.00% 

A12-1 

108 

109 

Closed 

No 

-33.70 

0.00 

-35.42 

4.86% 

A13-2 

108 

110 

Closed 

No 

-36.80 

0.00 

-39.58 

7.02% 

A14 

109 

111 

Closed 

Yes 

-27.80 

0.00 

-29.6 

6.08% 

A15 

109 

112 

Closed 

Yes 

-3.90 

0.00 

-4.35 

10.34% 

A16 

110 

111 

Closed 

Yes 

-19.80 

0.00 

-20.89 

5.22% 

A17 

110 

112 

Closed 

Yes 

3.90 

0.00 

4.35 

10.34% 

A18 

111 

113 

Open 

No 

0.00 

0.00 

0 

0.00% 

A19 

111 

114 

Closed 

No 

-56.60 

0.00 

-50.49 

10.80% 

A20 

112 

113 

Open 

No 

0.00 

0.00 

0 

0.00% 

A21 

112 

123 

Open 

No 

0.00 

0.00 

0 

0.00% 

A22 

113 

123 

Closed 

No 

-258.70 

7.81 

-265 

2.38% 

A23 

114 

116 

Closed 

No 

-50.40 

1.89 

-50.49 

0.18% 

A24 

115 

116 

Closed 

No 

204.90 

0.00 

204.49 

0.20% 

A25-1 

115 

121 

Open 

No 

0.00 

0.00 

0 

0.00% 

A26 

115 

124 

Closed 

No 

13.00 

0.00 

10.51 

19.15% 

A27 

116 

117 

Open 

No 

0.00 

7.59 

0 

0.00% 

A28 

116 

119 

Closed 

No 

313.50 

1.29 

309 

1 .46% 

A29 

117 

118 

Closed 

No 

122.20 

2.29 

123.55 

1.09% 

A30 

117 

122 

Closed 

No 

-122.20 

0.00 

-123.55 

1.09% 

A31-1 

118 

121 

Closed 

No 

-104.50 

0.54 

-104.73 

0.22% 

A32-1 

119 

120 

Closed 

No 

64.30 

0.01 

64 

0.47% 

A33-1 

120 

123 

Open 

No 

0.00 

0.08 

0 

0.00% 

A34 

121 

122 

Closed 

No 

-174.70 

0.00 

-176.45 

0.99% 

Table  2.  PWS’s  ACPF  versus  DCOPF:  IEEE  RTS  One  Area  case,  scenario  M=6. 

All  columns  are  as  explained  in  Table  1  with  the  exception  of  the  “Status” 
column  which  shows  which  lines  are  open.  Note  that  the  largest 
percentage  deviations  occur  across  lines  whose  flow  is  small. 

Results  for  scenarios  2  and  4,  summarized  in  table  3,  are  similar. 
Differences  in  line  power  flows  are  negligible  either  because  of  their  small 
relative  (percentage)  values,  or  because  of  the  low  power  flowing  across  lines. 
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Scenario 

Statistics  for  "Absolute  Deviation" 

Max.  (MW) 

Max.  (%) 

Average  (MW) 

Average  (%) 

Std.  dev. 
(MW) 

Std.  dev. 
{%) 

0 

8.04 

5.29 

1.70 

1.18 

1.96 

1.22 

2 

2.37 

13.04 

0.44 

0.92 

0.50 

2.28 

4 

10.37 

4.84 

1.24 

0.83 

2.10 

0.96 

6 

6.98 

19.15 

1.95 

3.68 

2.21 

0.96 

Table  3.  Overall  comparison  of  power  flows  provided  by  DCOPF  and  PWS.  The 
average  deviation  across  the  four  scenarios  is  less  than  5%. 

Overall,  the  differences  found  are  probably  acceptable  in  the  context  of 
solving  interdiction  problems.  This  concurs  with  the  generally  accepted  notion 
that  DCPF  renders  acceptable  results  in  the  context  of  contingency  analysis 
(e.g.,  Wood  and  Wollenberg  1996).  We  assume  that  the  differences  will  remain 
small  in  all  our  analyses. 
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IV.  BENDERS  DECOMPOSITION  FOR  THE  PROBLEM 
WITHOUT  SYSTEM  RESTORATION 


This  chapter  describes  Benders  decomposition  applied  to  the  MIP 
interdiction  model  LDI-DCOPF  without  system  restoration  presented  in  Chapter 
II,  Section  B.  The  chapter  begins  with  the  application  of  the  BDA  described  in 
Chapter  II,  Section  D  to  a  three-bus  test  case  with  interdictable  elements  limited 
to  lines.  The  general  case  is  then  addressed. 

A.  SMALL  TEST  CASE  PROBLEM  DERIVATION 
1.  Description 

This  section  illustrates  Benders  decomposition  applied  to  a  small,  three- 
bus  test  grid.  This  is  intended  to  familiarize  the  reader  with  the  mathematical 
derivation  for  the  general  case  that  will  be  described  later  in  this  chapter.  Even 
for  this  small  example,  the  mathematical  derivation  is  somewhat  tedious.  In 
order  to  simplify  the  exposition,  we  restrict  the  interdictable  components  to  lines 
only.  Figure  7  gives  the  one-line  diagram  for  the  example  grid. 


Bus  1  Bus  3 


Figure  7.  Three-bus  (and  three-line)  example.  Buses  1  and  2  are  connected  to  one 
generator  each.  Bus  3  is  connected  to  a  load  of  L.^  MW.  All  the  other 

symbols  in  the  figure  represent  decision  variables:  are  generator 

outputs;  6>,6>2,6»3  are  phase  angles  at  the  buses;  power 
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flows  on  the  lines. 


2.  Mathematical  Formulation 

In  this  example,  generating  costs  are  disregarded,  so  we  seek  to  minimize 
the  load  shed  at  bus  3.  The  problem  development  starting  at  the  DCOPF  level 
follows. 

Decision  variables(units): 

,  Generator  output  at  bus  i,  for  /  =  1,2  (MW) 

,  Power  flow  on  line  (/,  j) ,  for  (/,  y)  =  ( 1 ,2),  ( 1 ,3),  (2,3)  (MW) 

5*3 ,  Load  shed  at  bus  3  (MW) 

0^,  Phase  angle  at  bus  i,  for  i=  1,  2,  3  (radians) 

Problem  data  (units): 

L3 ,  Demand  at  bus  3  (MW) 

Line  resistance  and  reactance,  respectively,  for  (/,7)  =  (1,2),  (1,3),  (2,3)  (f2) 

T; 

By ,  Series  susceptance  of  line  (/,  j),  By  =  ^  ^ , 

ry  +Xy 

for(/,y)  =  (l,2),(l,3),(2,3)  (l/Q) 

P^ ,  Maximum  generating  capacity  at  bus  i  ,  for  i=  1,2  (MW) 

F’j',  Line  {i,j)  transmission  capacity,  for  {i,j)  =  (1,2),  (1,3),  (2,3)  (MW) 

(Recall  that  all  our  units  are  converted  to  p.u.  values  as  described  in 
Chapter  II,  Section  A). 

In  this  small  DCOPF  example,  we  seek  to  minimize  load  shed  at  bus  3, 

that  is: 
min 

Subject  to: 

Power  balance  constraints: 

P"-Pi-Pn  =0 

P2+Pn+P23  =0 
P,^,+P^,  +  S,  =L, 
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Admittance  constraints: 


Pi2-Bn^\  +  B^2^2  =0 

P\2~Pn^\  ■*■^3^3 

^23  “^23^2  +^23^3  =0 

Power  Generation  bound  constraints: 

P,^<Pf 

Line  capacity  constraints: 

pL  <  pL 
pL  y  _pL 
pL  <  pL 
pL  y  _pL 
pL  <  pL 

pL  ^  _pL 

^23  —  ''23 

Upper  bound  on  load-shedding  constraint; 

S,<L, 

Variable  sign  constraints: 

^3  >0 

P,^,  P^>Q 

Pi2,Pn,P22  unrestricted 
6^,  62, 6^  unrestricted 

We  can  now  extend  this  formulation  to  account  for  potential  interdiction  by 
introducing  three  variables  to  represent  interdiction  of  the  lines.  Note  that  the 
interdiction  variables  must  turn  off  or  turn  on  the  power  flow  on  the  lines  (i.e., 
open  or  close  the  lines,  respectively).  When  dy  =1,  the  line  {i,j)  is  interdicted, 

and  when  Sy  =  Q,  the  line  {i,j )  is  left  intact. 
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Letting  M.j=Py  +B.j6,j,  where  O.j  bounds  the  maximum  difference 
between  phase  angles  at  buses  i  and  j ,  the  interdiction  problem  for  given 
interdictions  variables  and  can  be  written  as  follows:  (Note:  The 

variables  tt  represent  the  dual  variables  of  the  balance  and  bound  constraints.) 
(I-DCOPF)  :  max  min  S. 

SeA 

s.t.: 


pG  pL 

M  M2  M3 

=  0 

(^f) 

P,^  +  P,'^,+P^, 

=  0 

(<) 

Pn  +  Pxi  +  *^3 

i^l) 

Pn  ~  P\2^\  P\2^2 

(<2) 

P\2  ~  P\2^\  P\2^2 

^-MnSn 

i^n) 

Pn  ~  Pn^\ 

<M„St, 

(<3) 

Pn  ~  P\2^\ 

(^if) 

^23  “  ^23^2  +  ^23^3 

<  M,Ji, 

(<3) 

P22  ~  ^23^2  ^23^3 

(<3) 

P^ 

<P° 

(<) 

P2 

<P° 

(<) 

pL 

12 

pL 

pL 

pL 

pL 

-*23 

</>‘(i-4)  (<■' 

^23 

*^3 

<L,  «”■') 

^3  >0 

>0 

Py2 ,  Py^ ,  .^23  unrestricted 

6^,6^,  6^  unrestricted 
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where: 


A  =  {-J, 


23 


<^12  +  <^13  +  ^ 


<  1 

23  —  -*-’^12 ’^13 ’^23  ' 


{0,1}) 


because  we  assume  one  unit  of  interdiction  resource  is  available. 

Notice  that  every  admittance  equation  has  been  split  into  two  inequalities. 
This  is  done  in  accordance  with  equation  (2.9)  (see  Chapter  II,  Section  B),  to 
linearize  the  multivariable  product  that  results  from  introducing  the  interdiction 
variables  into  the  model.  When  S^=0,  every  pair  of  inequalities  enforces 

+5.6>.  =  0.  When  dy  =1  (line  {i,j)  is  interdicted),  these  constraints  do 
not  bind  because  is  large.  However,  in  this  case  will  still  be  forced  to  be 
zero  by  constraints  If  and  If  >-^(\-Sl;) . 

The  optimal  interdiction  problem  consists  of  maximizing,  by  choice  of 
interdiction  variables  ,  the  amount  of  power  that  must  be  shed  in  the  system. 

As  shown  in  Chapter  II,  Section  B,  this  max-min  problem  can  be  easily 
reformulated  as  a  max-max  problem  by  taking  the  dual  of  the  inner  problem,  I- 
DCOPF,  listed  above.  This  yields: 

(DI-DCOPF):  max  max  + 

d  7t 

+3i,‘  (1  -  <y‘  )  +  />5  (1  -  )  +  L,  (2-3“ ) 

s.t. 

;rf+;rf  <0  (Z^'^) 

^2+<  ^0  {Pf) 

<1  (^3) 

-Ttl  +  7tl  +  <2  +  Ttf  +  =0  {Pff) 

-Ttf  +  7cl  -I-  ;rj^3  -I-  Ttf  -I-  -1-  =0  {Pf) 

+  ;r3"  +  <3  +  ;r2";  +  =  0  {P^, ) 
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~Px2^Xl  ~Px2^Xl  ~Px2^X2  ~ 

o 

II 

^Xl^Xl  ■*■  ^Xl^Xl  ~  -^23-^23  ~  ^12^22 

o 

II 

Px2^X2  Px2^X2  ^22^22  ^22^22 

o 

II 

^n^^X2^^22 

<1 

TT^ 

unrestricted 

<0 

„L*  „LCap* 

7t  ,  7t 

>0 

qL  qL  qL 

M2’  M3’  M3 

e  {0,1} 

The  resulting  model,  DI-DCOPF,  contains  a  number  of  cross-products  of 
the  form  5n  in  the  objective  function,  which  we  can  linearize  as  described  in 
Chapter  II,  Section  B.  This  linearization  yields  the  following  MIP: 

(LDI-DCOPF)  :  max  kI (v^  -  )  +  M^^  -  vf,  )  +  -  v^l ) 

oeA.v,^ 

+  ^2  )  +  ^3  ) 

■  pi  ( \—P^  A  _  A 

■'"''23  v^23  '^23  /  M2  M12  '^12  /  M3  M13  M3  / 

_pi/,,iCap'  _  iCap^\  ,  T  ^Load 
M3M23  M3  -^3'^3 


s.t. 

_iS  1  G 

<0 

^■5  , 

7^2  +7^2 

<0 

iP2) 

77^+71';°““ 

<1 

iS2) 

-Ttf  +77^  +7t^^  +7t^l  +77^2“'’*  +77^2“^' 

=  0 

iPn) 

-77I  +77^  +77^2  +^X2  +^n‘‘'’"  +77^2^’’ 

=  0 

iPx2) 

-77l  +77^  +<3  +77^2  +^22“”^  +^22“”' 

=  0 

(P22) 

—Bx277x2  —  B^2^^2  ”-^13-^13  ~  ^X2^X2 

=  0 

(^1) 

^X2^X2  ^X2^X2  ~  ^22^22  ~  ^22^22 

=  0 

(02) 

Px2^X2  Px2^X2  ^22^22  ^22^22 

=  0 

(03) 

^12 +  (^1^3 +  (^{3 

VG  V{d,77) 

<1 

where: 
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<0 


—L  —LCap  —G  L  LCap  —Load 

7t  ,  7t  ,  7t  ,  V  ,  V  ,  7t^ 

7t  ,  7t  ,  V  ,  V  >  (J 

TT^  unrestricted 


qL  cL 
^12  9  ^12  9  ^2 


G{0,1} 


The  linearization  of  constraints  ve  V(S,7r)  is  specified  in  detail  in  Appendix 


C,  but  as  an  example,  consider  the  constraint  =7r’:.  S^.  ,  where  >0.  The 

(/(/(/  y 


linearization  uses  the  following  constraints: 


V..  <  71:  O 

y  y  y 

<n,. 


_r  ^  =t 


^  ^0 


The  procedure  delineated  above  adheres  to  the  formulation  developed  by 
Salmeron  et  al.  (2004).  This  formulation  allows  us  to  use  any  general-purpose 
MIP  solution  technique  to  solve  LDI-DCOPF,  and  therefore,  l-DCOPF.  The 
remainder  of  this  section  is  devoted  to  demonstrating  the  application  of  Benders 
decomposition  to  LDI-DCOPF  in  our  three-bus  network,  following  the  steps  in 
Chapter  II,  Section  D. 


We  first  construct  the  subproblem  for  a  given  vector  of  complicating 
variables  4  iteration  k\ 


(SP^ {Sj^y) .  max .TTj  +  M^2 ^12  ) "*" (^13  ^13  ) '*' ^23 (^23  ^23 ) 

v,;T 

+  py^  )  +  ^3  -  ^13  ) 

4-P^  ( TT^^'^P  —  \  \  _  pi  A/Cap'  _  LCap+ \ 

■'■-'23  V‘'23  '‘'23  }  ^\1  v'^12  ^\1  J  '  13  V '^13  '^13  J 

_uP  f,,PCap'  _  LCap*  \  ,  j  Load 
''23  V  *^23  *^23  /■*■  ^3'^3 
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s.t. 

_L  -n-G 

71^  +  71^ 

< 

0 

(^4) 

—S  1  — G 

71  2  ■!■  ^2 

< 

0 

(^2^) 

,  ^Load 

71  2  -^71  2 

< 

1 

(*^3) 

+  <2 

+  <2 

.  ^LCap* 

-r  2*12 

.  —LCap' 

■*■2*12 

= 

0 

(^2) 

~7l^  -l-  71  2 

+  <3 

+  ^n 

,  —LCap^ 

+  71,2 

I  ^LCap 

T  2*13 

= 

0 

(^3) 

-7I2 

“*“•^23 

.  ^LCap* 

■'■2*23 

,  —LCap' 

T  2*23 

= 

0 

(^23) 

~^n^\2 

-B,271 

i-E 

'l3'^13  ~-®i3'^13 

= 

:0 

(^1) 

^\2^\2  ^\2^\2 

~  B22 

.'^23  ~  ^22^22 

= 

:0 

(^2) 

+ B22 

71  22  B 2^77 22 

= 

:0 

(^3) 

ve 

Constraints  ve  ^(4,-^)  are  linearized  as  in  Appendix  C  with  S  replaced 

by  4- 


The  master  problem  at  iteration  k  becomes: 


(MPj^):  max  z 


<yeA,ZeM 

s.t 


Z<  ^i2^\2Yl2l,k'  '^12  ^12)?n3,  ■^13  ^13?ni,/t'  '^13  ^13)?n3,i:'  ■^23  ^23?^31,t' 

—^23  ~  ^2'i)y2'i'i,k'  ~^\1  ^\lY\l\,k'  '^12  ~  ^12)^123,*'  “'^13  ^13^131,*’  '^13  ~  ^13)?n3,i:' 

—TT^S  (l  —  S  ^Y^  S  _jfLCap  \ykCap  x  jfLCap  s  yLCap 

'*'23  ^23/  231, t'  ^ '*'23  ^23//  233,*’  ^ '*'12  ^12/121,*'  '*'12  ^122/123,*’  ^'*'13  m3/ 131,*' 

_~LCap*  i-x  _  ^  \  LCap*  ,  =LCap* ^  ^.LCap*  _  =LCap* ^  ^^J-Cap*  _  =LCap^ ^  ..LCap^ 

2*13  ^132/133,*'  "''^^23  ^23/  231,*'  ^^23  ^232/  233,*'  ^^12  ^12/121,*’ 

,=LCap^  /IP  \^Cap^  _  =LCap^  ^  ^Cap^  ,  =LCap^  *1  —  /?  ^  ^Cap^ 

■^2*12  U  ^122/123,*'  2*13  ^13/131,*'  “^2*13  1,1  ^132/133,*'  2*23  ^23/  231,*' 

+•^23  ^  (1  ~  ^23)?M3,*'  “*■ ‘^3,*’  “*■  •^12  ^12,*'  “*■  •^13  ^13,1  “*■  •^23  ^23,*'  ~’^12  ^12,*'  ~’^13  ^13,*'  “'^23^23,*' 


,=LCap*  LCap*  ,  =LCap*  LCap*  ,  =LCap*  LCap*  _  =LCap  „LCap  _  =LCap  „LCap  _  =LCap  „LCap 

-r2*i2  '/12,*'  ■'■2*13  '/l3,*'  ■'■2*23  '/23,*'  2*12  'lu,k'  2*13  '/13*'  2*33  I/2: 


/23,*' 


VA:'  =  {1,2,...,A:} 

In  the  above  formulation,  the  subscript  ^for  any  variable  {e.g.,Yl2ik')  refers 


to  the  optimal  value  of  the  incumbent  variable  provided  by  the  solution  to 
(SP^,(4’))  at  iteration  k'\k'<k. 


Having  developed  our  master  problem  and  subproblem,  we  implement  the 
BDA  in  Chapter  II,  Section  D,  and  test  the  three-bus  case  for  convergence.  We 
use  the  following  values  for  the  test  case: 
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L^=5  p.u. 

5^.  =5p.u. 

^^P/  =  3p.u. 

^2, ^3, ^23  =4  p.u. 

The  algorithm  converges  successfully  to  the  same  optimal  interdiction 
plan;  (42  =  0,43  =  1^43=0)  that  of  the  MIP  formulation  in  LDI-DCOPF.  The 
following  graph  illustrates  how  upper  and  lower  bounds  converge  after  four 
iterations  of  the  algorithm. 


BOUNDS  VERSUS  ITERATION 


Figure  8.  Convergence  of  Benders  decomposition  for  the  three-bus  case.  The 
optimal  interdiction  plan,  which  consists  of  interdicting  line  (1,3)  in  Figure 

7,  sheds  a  load  of  one  p.u.. 

This  section  has  illustrated  Benders  decomposition  using  a  small  example. 
This  foundation  should  help  the  reader  understand  the  generalization  of  the 
procedure  to  any  power  grid  in  which  buses,  generators,  substations  and  lines 
can  be  interdicted.  The  extension  to  incorporate  the  effect  of  system  restoration 
over  time  will  be  addressed  in  Chapter  V. 
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B.  BENDERS  DECOMPOSITION:  GENERAL  INTERDICTION  PROBLEM 

WITHOUT  SYSTEM  RESTORATION 

The  general  interdiction  problem  should  be  scalable  to  let  us  analyze  an 
arbitrary  power  system.  We  will  not  discuss  the  derivation  of  the  BDA  for  I- 
DCOPF  as  done  in  the  previous  section,  however.  Appendix  A  provides  a 
complete  derivation  of  the  Benders  decomposition  procedure  for  l-DCOPF.  The 
derivation  is  broken  down  into  these  steps: 

A.1  shows  the  interdiction  model,  l-DCOPF. 

A.2  shows  the  dual  of  the  interdiction  model,  DI-DCOPF. 

A.3.  shows  the  dual  of  the  interdiction  model  after  linearization  of  5n 
products,  LDI-DCOPF. 

A.4  shows  the  BDA  master  problem. 

The  subproblem  at  iteration  k  is  LDI-DCOPF  with  all  ^’s  replaced  by 
fixed  4  ’s  (given  by  the  user  at  iteration  ^  =  0  and  by  {MPk)  at  iteration  k  >  0), 
where  4  g  A,VA: . 

Once  the  subproblem  is  solved  for  a  given  interdiction  plan,  we  retrieve 
the  dual  values  for  each  subproblem  constraint  and  use  them  as  coefficients  for 
the  master  problem.  The  master  problem  at  iteration  k  is  displayed  in  Appendix 
A.4. 

C.  RESULTS 

The  BDA  is  implemented  in  GAMS  (GAMS  2003)  and  solved  with  CPLEX 
version  8.1  (GAMS-CPLEX  2003)  as  the  underlying  solver,  and  run  on  a  Pentium 
4  laptop  computer  at  3.0  GHz  and  with  512  Mbytes  RAM.  We  turn  off  the 
“presolve”  option  in  CPLEX  because  it  leads  us  to  erroneous  duals  from  the 
subproblem. 

Figure  9  shows  the  sequence  of  lower  and  upper  bounds  generated  by 
BDA  for  the  IEEE  RTS  One  Area  case  (1999-1),  with  interdiction  resource  M  =  6 
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and  the  assumptions  made  in  Salmeron  et  al.  (2003).  The  optimal  interdiction 
plan  for  this  case  matches  the  solution  obtained  by  solving  the  MIP  LDI-DCOPF 
directly. 


BOUNDS  VERSUS  ITERATION 


Figure  9.  IEEE  RTS  One  Area  case  without  system  restoration:  Convergence.  Note 
the  significant  improvement  of  the  lower  and  upper  bound  at  early 
iterations.  As  the  number  of  iterations  increases  the  effectiveness  of  the 
cuts  (given  by  the  upper  bound  from  the  master  problem)  decreases. 
Initial  upper  bound  is  140,000,  but  the  vertical  axis  on  the  graph  is 
truncated  for  display  purposes. 

The  results  for  this  problem  show  that,  by  iteration  375,  the  lower  bound  is 
already  within  7.6%  of  the  optimal  solution,  while  the  upper  bound  is  within  12%. 
(For  this  comparison  we  use  the  optimal  solution  obtained  with  LDI-DCOPF.)  We 
achieve  these  bounds  in  approximately  10  minutes;  however,  it  takes  1  hour  and 
20  minutes  to  reach  the  optimal  solution  which  is  14,048  $/hr.  This  result  shows 
Benders  decomposition  may  obtain  sensible  bounds  in  a  relatively  small  number 
of  iterations,  and  an  acceptable  time.  However,  reaching  optimality  takes 
considerable  effort.  This  difficulty  will  be  addressed  in  Chapter  VI. 
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Figure  10  shows  how  the  time  required  to  solve  the  master  problem  and 
subproblem  change  as  the  algorithm  proceeds.  This  figure  reveals  the 
increasing  difficulty  of  the  master  problem  as  more  cuts  are  added  with 
subsequent  iterations.  While  the  time  to  solve  each  succeeding  master  problem 
of  the  decomposition  process  significantly  increases,  the  relative  efficiency  of  the 
new  cuts  are  less  significant  as  shown  by  Figure  9. 
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Figure  10.  IEEE  RTS  One  Area  case  without  system  restoration:  Time  versus 
Iteration  .  The  time  to  solve  each  master  problem  in  the  algorithm 
significantly  increases  with  subsequent  iteration.  The  maximum  time  to 
solve  any  subproblem  is  less  than  0.2  seconds  while  master-problem 
solution  times  are  as  high  as  18  seconds. 
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V.  BENDERS  DECOMPOSITION  FOR  THE  PROBLEM  WITH 

SYSTEM  RESTORATION 


This  chapter  extends  Benders  decomposition  to  the  interdiction  problem 
with  system  restoration.  We  describe  the  models  involved  in  the  new 
decomposition,  and  show  preliminary  results.  (Chapter  VI  will  show  improved 
results  through  some  refinements  in  BDA.) 

A.  FORMULATION 

This  section  illustrates  the  application  of  Benders  decomposition  to  the 
interdiction  problem  with  system  restoration,  l-DCOPF-R.  The  derivation  of  I- 
DCOPF-R  is  included  in  Appendices  B.1  and  B.2.  We  will  not  discuss  the 
complete  derivation  of  the  Benders  decomposition  procedure  for  the  system 
restoration  case,  which  is  included  in  Appendix  B.  Specifically: 

B.1  shows  time-period  constructs. 

B.2  shows  the  interdiction  model  l-DCOPF-R. 

B.3  shows  DI-DCOPF-R,  which  is  the  interdiction  model  with  the  inner 
power  flow  model  replaced  by  its  dual. 

B.4  shows  DI-DCOPF-R  after  linearization  of  cross-products.  We  denote 
this  model  LDI-DCOPF-R. 

Finally,  B.5  shows  the  BDA  master  problem. 

We  start  the  description  after  the  linearization  of  the  dn  cross-products. 
That  is,  assuming  we  have  LDI-DCOPF-R  as  the  MIP  to  which  we  will  apply 
Benders  decomposition. 

Again,  the  subproblem  at  iteration  k  is  the  LDI-DCOPF-R  with  all  5's 
replaced  by  fixed  4  ’s  (given  by  the  user  at  iteration  k  =  0  and  by  the  {MPk)  at 
iteration  k>0),  where  Si^eA,\/k,  and  the  dual  values  for  each  subproblem 
constraint  are  used  as  coefficients  for  the  master  problem. 


45 


The  BDA  master  problem  is  shown  in  Appendix  B.5.  At  each  iteration  k,  a 
new  cut  is  added  to  the  master  problem.  This  provides  a  (monotonically 
decreasing)  upper  bound  z^,  and  an  interdiction  plan  that  will  be  used  in  the 

subproblem  (SP^(4))  to  calculate  a  lower  bound. 

B.  RESULTS 

Figure  11  shows  the  convergence  of  BDA  for  the  IEEE  RTS  One  Area 
case  with  system  restoration.  Note  the  significant  improvement  in  the  bounds 
during  the  early  iterations  of  the  decomposition  process  when  compared  to  the 
last  iterations. 


BOUNDS  VERSUS  ITERATION 


Iteration 


Figure  1 1 .  IEEE  RTS  One  Area  case  with  system  restoration:  Convergence.  The 
upper  and  lower  bounds  of  this  problem  converge  quickly  in  early 
iterations.  Convergence  is  reached  in  120  iterations  and  takes  33 
seconds.  Also  note  that  the  optimal  objective  value  (circled  in  graph)  is 
reached  in  early  iterations,  but  it  takes  the  algorithm  a  relatively  long  time 
to  prove  it  has  obtained  an  optimal  solution. 
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Although  the  case  with  system  restoration  is  intuitively  more  difficult  to 
solve  than  the  case  without  system  restoration  (because  of  the  number  of  time 
periods,  and  the  related  number  of  additional  variables),  the  former  actually 
solves  faster  in  the  IEEE  RTS  One  Area  case.  This  could  be  attributed  to  the 
fact  that,  for  the  problem  with  restoration,  some  candidate  components,  such  as 
those  with  the  largest  repair  times,  are  more  obviously  attractive  than  others  to 
become  part  of  the  optimal  solution.  This  might  allow  the  BDA  to  target  solutions 
that  include  those  components  at  early  stages,  which  in  turn  provides  accurate 
bounds  sooner  than  in  the  case  without  system  restoration  (in  which  the 
attractiveness  of  all  components  is  more  balanced).  Our  results  support  this 
conjecture.  While  it  takes  over  700  iterations  and  1  hour  and  20  minutes  to  solve 
the  IEEE  RTS  One  Area  case  without  system  restoration  problem,  the  problem 
with  system  restoration  solves  in  only  120  iterations  and  33  seconds.  We  find 
similar  behavior  in  the  RTS  Two  Area  case. 

Figure  12  shows  that  master  problem  time  by  iteration  remains  almost  the 
same  for  the  IEEE  RTS  One  Area  case  with  system  restoration.  We  will  show 
that  this  is  not  the  case  when  we  apply  the  BDA  to  the  IEEE  RTS  Two  Area  case. 


TIME  BY  ITERATION 


Figure  12.  IEEE  RTS  One  Area  case  with  system  restoration:  Cumulative  time.  Time 
is  shown  for  the  subproblem  and  for  the  master  problem  separately.  Note 
that  the  time  to  solve  the  master  problem  appears  to  remain  stable  by 
iteration,  despite  the  increasing  number  of  constraints. 
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To  illustrate  the  potential  difficulty  in  solving  the  master  problem  as  the 
number  of  iterations  increases,  we  analyze  the  results  for  interdiction  scenario 
M  =  12  for  the  IEEE  RTS  Two  Area  case.  Figure  13  shows  the  results  for  a 


2,000-iteration  run  of  this  problem. 
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Figure  13.  IEEE  RTS  Two  Area  case  with  system  restoration:  Solution  trajectory.  In 
2,000  iterations,  taking  75  hours,  a  gap  of  16%  is  achieved.  The  top  left 
graph  shows  that  a  near-optimal  solution  is  found  at  the  early  stages,  but  it 
takes  many  iterations  to  prove  it.  The  top  right  graph  shows  the 
cumulative  solution  time,  by  iteration,  for  the  master  problem  and 
subproblem.  Note  how  the  master-problem  solution  times  increase 
substantially  as  the  algorithm  proceeds.  The  bottom  left  plot  shows  how 
the  optimality  gap  changes  with  time.  Note  the  large  decrease  in  the  gap 
in  the  early  iterations.  The  bottom  right  picture  is  an  expanded  view  of  the 
first  30  minutes  of  the  “GAP  VERSUS  TIME”  plot.  The  gap  is  reduced  to 
50%  in  the  first  30  minutes.  The  problem’s  lower  bound  is,  in  actuality, 
only  2%  from  the  optimal  solution  value  at  this  time,  although  the  algorithm 
cannot  prove  this  within  the  first  75  hours. 


Clearly,  BDA  is  impractically  slow  for  this  problem.  This  suggests  the 
need  for  strategies  to  accelerate  BDA’s  convergence.  The  next  chapter 
describes  and  demonstrates  techniques  we  have  developed  to  do  this. 
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VI.  ALGORITHM  REFINEMENTS  AND  RESULTS 


As  seen  in  Chapter  V,  long  master-problem  solution  times  are  a  major 
factor  in  the  overall  efficiency  (or  inefficiency)  of  BDA.  This  prompts  us  to  focus 
on  the  master  problem  and  reduce  its  solution  times  to  accelerate  overall 
convergence  of  BDA. 

A.  INTRODUCTION 

This  chapter  explores  the  following  techniques  to  reduced  master-problem 
solution  times: 

1.  Solving  a  relaxed  master  problem  in  some  iterations  rather  than  the 
standard  master  problem, 

2.  Dropping  certain  Benders  cuts  as  iterations  proceed,  and 

3.  Not  solving  the  master  problem  to  optimality  in  all  iterations  (Sub- 
optimal  integer  solutions). 

The  modified  BDA  which  includes  the  above  strategies  1,  2  and  3  is 
shown  in  Figure  14. 

The  various  techniques  and  associated  results  are  discussed  in  Sections 
B  through  D  below.  A  combined  strategy  is  reviewed  in  Section  E. 
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Figure  14.  Benders  Decomposition  Algorithm  (with  refinements)  flowchart.  This 
diagram  depicts  the  steps  of  our  enhanced  BDA,  incorporating  the 
proposed  strategies  to  speed  up  convergence.  RMP  refers  to  the  linear 
relaxation  of  the  BDA  master  problem.  The  flowchart  includes  the 
refinement  techniques  that  are  explained  in  Sections  B  through  D  of  this 

Chapter. 
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B.  RELAXING  THE  MASTER  PROBLEM 

The  first  technique  we  try  solves  a  relaxed  master  problem  (RMP^)  in  most 
iterations  k,  and  only  solves  the  true  (mixed-integer)  master  problem  periodically. 
A  relaxed  master  problem  should  be  much  faster  to  solve,  and  its  use  may  lead 
to  an  overall  reduction  in  solution  time.  We  refer  to  this  technique  as  the 
“relaxed-MP  strategy.” 

RMPyt  is  formed  by  treating  any  discrete  variables  as  continuous  (while  the 
rest  of  the  problem  remains  the  same).  The  optimal  value  of  a  relaxed 
maximizing  model  yields  an  upper  bound  on  the  optimal  value  of  the  full  model. 
However,  the  solution  to  the  relaxed  master  problem  is  not  a  feasible  interdiction 
plan  (unless  it  happens  to  be  an  integer  solution).  For  this  reason,  the  relaxed 
master  problem  solution  cannot  be  used  in  the  subproblem  to  obtain  a  valid  lower 
bound;  however,  the  Benders  cut  generated  by  the  subproblem  at  that  solution  is 
valid.  In  order  to  improve  the  lower  bound,  we  must  solve  the  actual  master 
problem,  and  we  do  this  at  a  specified  interval  (for  example,  every  ten  iterations). 
Unless  specified  otherwise,  in  the  test  cases  that  follow,  our  relaxed-MP  strategy 
consists  of  solving  the  true  master  problem  once  every  ten  iterations  and 
otherwise  solving  RMPyt. 

To  show  the  benefits  of  using  the  relaxed-MP  strategy,  we  first  consider 
BDA  applied  to  the  IEEE  RTS  One  Area  case  without  system  restoration.  The 
original  algorithm’s  solution  trajectory  for  this  problem  is  presented  graphically  in 
Chapter  IV,  Figure  9.  Recall  that  it  originally  takes  about  1  hour  and  20  minutes 
to  solve  this  case.  The  relaxed-MP  strategy  solves  the  problem  in  only  26 
minutes,  a  67%  improvement  over  the  original  algorithm.  Figure  15  shows  how 
the  problem  converges. 

The  reduction  in  overall  solution  time  cannot  be  attributed  solely  to  faster 
solutions  of  the  master  problem,  however,  since  the  modified  algorithm  solves  in 
only  281  iterations  compared  to  the  original  768.  It  appears  that  cuts  generated 
from  non-integer  solutions  may  be  more  effective  than  standard  cuts  generated 
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at  integer  solutions.  We  are  unsure  why  this  occurs  in  this  case,  and  in  fact,  we 
will  see  later  that  this  behavior  cannot  be  generalized  to  all  our  cases. 


Figure  15.  Relaxed  master  problem  strategy:  Convergence  (I).  IEEE  RTS  One  Area 
case  without  system  restoration.  We  solve  the  full  master  problem  every 
10**^  iteration  and  otherwise  solve  that  problem’s  continuous  relaxation. 
This  problem  solves  67%  faster  using  this  technique  compared  to  the 
original  BDA  which  solves  the  (mixed-integer)  master  problem  in  each 

iteration. 

Figure  16  shows  the  solution  times  by  iteration,  for  the  master  problem 
and  for  the  subproblem. 
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Figure  16. 
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Relaxed  master  problem  strategy:  Time  versus  iteration  (I).  IEEE  RTS 
One  Area  case  without  system  restoration.  We  depict  computational 
times  when  using  the  relaxed-MP  strategy.  The  top  figure  shows  the 
subproblem  solution  time  by  iteration  (which  averages  under  0.15 
seconds).  The  full  master  problem  is  solved  every  lO*'^  iteration  (these 
are  easily  distinguished  in  the  bottom  plot  due  to  their  longer  solution 
times).  The  average  time  required  to  complete  each  master  problem 
iteration  is  reduced  to  less  than  one  second  in  contrast  to  the  original  BDA 
whose  iteration  average  is  9.5  seconds. 


Figure  16  shows  that  subproblem  solution  times  are  approximately  the 
same  as  in  the  original  BDA  (see  Chapter  IV,  Figure  10),  whereas  the  master 
problem  solution  times  have  been  reduced  considerably.  The  apparently 
exponential  increase  in  master-problem  solution  times  (exhibited  by  the  original 
BDA  as  iterations  progressed)  does  not  appear  now:  The  new  master-problem 
solution  times  appear  to  increase  only  linearly  by  iteration.  The  average  time  per 
iteration  using  the  relaxed-MP  strategy  is  0.15  seconds,  much  less  than  the 
average  of  9.5  seconds  per  iteration  for  the  original  BDA. 


The  results  above  are  for  the  problem  without  system  restoration.  Figure 
17  depicts  the  solution  trajectory  for  the  IEEE  RTS  One  Area  case  with  system 
restoration,  using  the  relaxed-MP  strategy. 
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Figure  17.  Relaxed  master  problem  strategy:  Convergence  (II).  IEEE  RTS  One  Area 
case  with  system  restoration.  The  algorithm  closes  the  gap  rapidly  in 
early  iterations  but  takes  many  iterations  to  prove  optimality.  Detailed 
solving  time  from  the  algorithm,  not  shown,  reveals  that  a  10%  gap  is 
reached  in  12  seconds,  but  optimality  is  not  reached  until  60  seconds. 

Note,  however,  that  this  problem  takes  longer  to  reach  optimality  using  the 
relaxed-MP  strategy  than  the  original  BDA.  This  fact  can  be  attributed  to  the 
number  of  relaxed  master  problems  (nine  in  this  run)  that  are  solved  before  we 
solve  the  standard  master  problem  in  order  to  update  the  lower  bound.  If  we 
reduce  the  interval  for  solving  the  standard  master  problem,  for  example,  every 
3'^'^  iteration,  we  would  be  able  to  solve  this  problem  to  optimality  in  only  23 
seconds,  a  reduction  of  8  seconds  (25%)  over  the  original  BDA. 

Figure  18  shows  the  cumulative  solution  time  for  the  master  problem  for 
both  the  original  BDA  and  the  relaxed-MP  strategy.  Note  the  reduced  solution 
time  for  the  master  problem  when  compared  to  the  original  BDA. 


54 


Figure  18.  Relaxed  master  problem  strategy:  Time  versus  iteration  (II).  IEEE  RTS 
One  Area  case  with  system  restoration.  We  show  cumulative  time  to 
solve  the  master  problem.  The  average  time  per  iteration  for  the  master 
problem  in  the  relaxed-MP  algorithm  is  0.03  seconds  (including  iterations 
where  the  master  problem  is  solved  exactly),  a  substantial  improvement 
from  the  original  BDA,  which  is  0.15  seconds.  Cumulative  time  is  reduced 

accordingly. 

The  advantages  of  the  relaxed-MP  algorithm  are  more  significant  as 
problems  become  more  dificult.  This  algorithm  helps  close  the  optimality  gap 
faster  because  (1)  every  relaxed  Benders  cut  is  valid  and  can  therefore 
contribute  to  the  reduction  of  the  upper  bound,  and  (2)  the  relaxed  master 
problem  is  so  much  faster  to  solve.  Figure  19  shows  the  effect  of  this  algorithmic 
strategy  for  the  IEEE  RTS  Two  Area  case  with  system  restoration. 
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Figure  19.  Relaxed  master  problem  strategy:  Convergence  and  Time  versus 

Iteration  (III).  IEEE  RTS  Two  Area  case  with  system  restoration.  The  top 
graph  illustrates  that  the  relaxed-MP  algorithm  closes  the  optimality  gap 
more  in  earlier  iterations.  The  bottom  graph  illustrates  the  master 
problem’s  reduced  cumulative  solution  time  when  we  use  the  relaxed-MP 
algorithm.  Note  that  a  logarithmic  scale  on  the  time  axis  is  used  for  clarity. 

In  the  IEEE  RTS  Two  Area  case  we  can  appreciate  the  impact  of  the 
relaxed-MP  strategy  on  overall  solution  time.  The  original  BDA  takes  over  75 
hours  to  complete  2,000  iterations,  reaching  a  16%  gap.  By  solving  the  relaxed 
master  problem  for  nine  of  ten  iterations,  and  solving  the  standard  master 
problem  every  10**^  iteration,  we  reduce  the  total  solution  time  to  a  mere  26 
minutes  for  2,000  iterations,  and  reach  a  19%  gap.  The  original  algorithm 
reaches  a  20%  gap  at  iteration  1,683  in  approximately  37  hours;  in  contrast,  the 
relaxed-MP  algorithm  achieves  the  same  gap  in  1,762  iterations  but  only  in  20 
minutes.  The  average  time  per  iteration  for  this  test  case  using  the  original 


56 


algorithm  is  two  minutes  and  15  seconds;  using  the  relaxed-MP  algorithm,  we 
reduce  the  average  time  per  iteration  to  0.78  seconds.  The  original  BDA,  of 
course,  takes  much  more  time  to  achieve  such  gap. 


C.  CUT-DROPPING  STRATEGIES 

The  relaxed-MP  algorithm  helps  solve  certain  problems,  but  larger 
problems  need  additional  techniques  like  “cut-dropping”  if  they  are  to  be  solved 
efficiently.  The  goal  of  cut-dropping  is  to  limit  the  number  of  Benders  cuts  in  the 
problem  so  that  the  master  problem  remains  efficiently  solvable.  Of  course, 
convergence  of  BDA  may  be  lost  unless  care  is  taken. 

We  explore  the  following  cut-dropping  techniques: 

-  Drop  (delete)  all  “sufficiently  slack  cuts,” 

-  Drop  all  “non-LP-active  cuts,” 

-  Dropping  the  first  slack  cut, 

-  Keep  a  minimum  number  of  cuts  plus  all  active  cuts,  and 

-  Keep  the  n-most  active  cuts. 

The  following  subsections  explain  and  demonstrate  each  technique  in 

detail. 

1.  Dropping  All  Non-active  Cuts 

This  strategy  attempts  to  keep  only  binding  cuts  (the  cuts  that  remain 
“active”  between  iterations)  in  the  master  problem.  In  a  linear  problem,  an 
inequality  constraint  is  deemed  “active”  (or  binding)  at  a  given  feasible  solution  if 
the  exact  equality  holds  at  that  solution.  A  non-active  constraint  is  said  to  be 
slack.  The  concept  of  slackness  is  important  here  because  though  basic 
sensitivity  analysis  one  can  show  that,  at  the  optimal  solution  to  a  linear  problem, 
slack  constraints  can  be  dropped  without  changing  the  optimal  solution  to  the 
problem.  Moreover,  it  is  well  known  that,  for  linear  problems  solved  through 
Benders  decomposition,  one  may  eliminate  non-active  cuts  in  the  master 
problem  at  a  given  iteration  without  losing  convergence  to  an  optimal  solution. 
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For  a  mixed-integer  master  problem,  the  above  property  is  also  true,  but 
checking  for  active  cuts  is  a  more  complicated  task  because  a  cut  might  be  (in 
fact,  will  be  in  most  cases)  slack  at  the  optimal  integer  solution,  yet  removing  it 
could  cause  a  (strict)  relaxation  of  the  master  problem.  Checking  for  active  cuts 
in  a  MIP  requires  removing  the  cut  and  solving  the  new  MIP  (in  order  to  assess  if 
the  optimal  objective  function  value  has  changed).  Doing  this  for  every  cut  is  not 
likely  to  yield  an  efficient  algorithm,  of  course.  Instead,  our  approach  consists  of 
dropping  cuts  that  are  “estimated”  to  be  non-active.  We  accomplish  this  by 
removing  all  cuts  at  any  iteration  that  have  a  slack  greater  than  a  percentage  of 
the  optimal  value  of  the  current  MP.  The  slack  for  a  cut  of  the  form  z<c  -5 
(where  zgM,  c,^eM")  at  a  given  feasible  solution  {z,d)  is  given  by 
s  =  c*  ■d-z>Q.  Then,  5  is  compared  with  the  master-problem  objective  function 
value  z.  Ideally,  we  would  use  {z,d)  =  {z  ,d*)  (optimal  solution  to  the  master 
problem)  to  evaluate  the  (relative)  cut  slackness.  However,  sometimes  we  must 
rely  on  relaxed-MP  solutions  (see  Section  B)  or  even  on  sub-optimal  solutions 
(see  Section  D)  to  carry  out  this  comparison. 

We  use  the  RTS  One  Area  case  for  basic  tests  and  perform  several  runs 
with  various  “slackness  criteria.”  In  particular,  we  eliminate  any  cut  whose  slack 
s  is  greater  than  0.1%,  1%,  5%,  10%  and  20%,  respectively  of  the  current  z. 
The  results  for  all  these  cases  are  similar  and  unsatisfactory:  It  appears  that  this 
technique  eliminates  some  cuts  that  are  necessary  for  the  convergence  of  the 
problem.  For  example,  if  we  use  5>0.20z  (20%),  the  master  problem  bound 
does  not  improve  after  the  85*^^  iteration,  impeding  this  variant  of  BDA  from 
converging  (see  Figure  20). 
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Figure  20.  Dropping  non-active  cuts:  Convergence.  IEEE  RTS  Two  Area  case 
without  system  restoration.  Lower  and  Upper  Bounds  versus  number  of 
iterations.  We  only  use  estimated  active  cuts  in  the  master  problem. 

These  cuts  cannot  close  the  optimality  gap. 

2.  Dropping  Non-LP  Active  Cuts 

This  strategy  estimates  the  active  cuts  in  the  master  problem  by  using  the 
active  cuts  in  the  relaxed  master  problem.  We  keep  only  LP-active  cuts  in  the 
master  problem.  We  consider  a  cut  to  be  LP-active  if  it  is  active  in  the  relaxed 
master  problem.  This  technique  showed  no  significant  improvement  over  the 
previous  cut  dropping  technique  because  it  also  appears  to  remove  too  many 
cuts  from  the  master  problem.  For  example,  in  the  IEEE  RTS  One  Area  case 
with  system  restoration,  only  four  LP-active  cuts  remain  in  iteration  566,  yielding 
a  weak  upper  bound. 

3.  Dropping  the  First  Slack  Cut 

We  check  all  incumbent  cuts  at  any  given  iteration  until  we  find  one  whose 
slack  exceeds  a  pre-specified  threshold,  for  example,  being  20%  above  the 
master  problem  objective  function  value.  If  none  of  the  cuts  satisfies  this 
criterion,  none  is  dropped.  In  addition,  a  cut  that  is  removed  is  no  longer 
considered  in  the  master  problem. 
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This  strategy  still  removes  cuts  that  are  needed  for  convergence.  In  the 
IEEE  RTS  One  Area  case  with  system  restoration  and  a  slackness  threshold  of 
20%,  only  nine  cuts  remain  in  the  master  problem  by  iteration  300  and  the 
optimality  gap  is  597%.  When  we  increase  the  criterion  to  80%,  52  cuts  remain 
at  iteration  300,  but  the  optimality  gap  is  still  unacceptable  at  180%. 

A  variant  of  this  strategy  consists  of  selecting  the  first  cut  whose  dual 
variable  (at  the  optimal  branch-and-bound  node  for  the  master  problem  solution) 
is  zero  for  elimination.  Again,  this  alternative  does  not  provide  us  with  the 
necessary  set  of  cuts  to  close  the  optimality  gap. 

4.  Keeping  the  n>Most  Active  Cuts 

The  results  in  the  previous  subsections  show  that  we  must  sustain  an 
elevated  number  of  cuts  in  order  to  avoid  destroying  convergence.  Additionally, 
we  can  see  from  those  strategies  that  if  we  eliminate  cuts  based  on  a  pre-set 
slackness  criterion  we  may  be  eliminating  cuts  needed  for  convergence. 

Here  we  try  a  new  strategy  where  we  limit  the  number  of  cuts  to  a  pre¬ 
specified  value,  n.  We  assume  the  master  problem  with  n  cuts  is  solvable  in  a 
reasonable  amount  of  time.  At  every  iteration,  we  replace  the  cut  with  largest 
slack  by  the  cut  generated  at  the  incumbent  iteration.  This  strategy  guarantees 
that  the  “best”  n-^  cuts  are  always  used,  keeping  the  problem  to  a  manageable 
size  that  will  solve  relatively  quickly. 

Testing  indicates  that  being  too  restrictive  in  the  value  given  to  n  can 
prevent  the  algorithm  from  converging.  For  example.  Figure  21  shows  that  for 
the  IEEE  RTS  One  Area  case  with  system  restoration,  restricting  the  number  of 
cuts  to  n  =  30  is  not  efficient,  whereas  the  algorithm  converges  successfully  for  n 
=  110. 
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Figure  21 .  Keeping  the  n  most  active  cuts:  Convergence.  IEEE  RTS  One  Area  case 
with  system  restoration.  Gap  versus  iteration  plot.  Limiting  the  number  of 
cuts  makes  the  problem  simpler  to  solve.  However,  as  illustrated  in  this 
figure,  too  few  cuts  will  prevent  the  problem  from  closing  the  optimality 

gap. 

Specifically  for  n  =  30  cuts,  we  can  only  close  the  gap  to  60%  after  180 
iterations  (in  45  seconds),  and  cannot  improve  after  then.  For  n  =  60  cuts,  we 
can  achieve  a  34%  gap  in  62  iterations  (19  seconds).  For  n  =  90  cuts,  a  gap  of 
4%  is  achieved  in  99  iterations  (30  seconds).  Finally,  for  any  n  =  110  cuts  or 
more,  the  problem  solves  to  optimality  in  117  iterations  (35  seconds). 


D.  SUB-OPTIMAL  INTEGER  SOLUTIONS 

In  this  section  we  consider  a  strategy  that  is  similar  to  the  master-problem 
relaxation  described  in  Section  B.  This  extension  involves  establishing  a 
termination  criterion  for  the  (full,  mixed-integer)  master  problem  before  it  is 
solved  to  optimality.  Typical  termination  criteria  are  based  on  limiting  any  of  the 
following:  the  number  of  integer  solutions  found  during  the  branch-and-bound 
process,  the  number  of  nodes  explored  in  the  branch-and-bound  tree,  the 
computational  time  spent,  or  combinations  of  the  above,  such  as:  “Stop  after  a 
maximum  number  of  seconds,  if  at  least  one  integer  solution  has  been  found.” 
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If  the  integer  solution  found  is  not  integer-optimal  (due  to  the  early 
termination  according  to  pre-specified  criteria)  then  the  upper  bound  cannot  be 
updated.  However,  it  is  still  a  candidate  solution  for  improving  the  lower  bound, 
which  will  be  obtained  after  solving  the  subproblem  for  that  candidate  solution. 

To  ensure  that  the  upper  bound  eventually  improves,  at  least  every  m 
iterations  (e.g.,  m  =  50)  the  full  master  problem  needs  to  be  solved  to  optimality. 
(Recall  that  the  upper  bound  can  also  be  updated  whenever  we  solve  a  relaxed 
master  problem).  Figure  22,  shows  the  bound  trajectory  graph  for  the  RTS  Two 
Area  case  where  (a)  the  master  problem  is  solved  to  optimality  every  50**^ 
iteration,  (b)  the  master  problem  is  solved  to  a  sub-optimal  integer  solution  every 
10*^  iteration  that  is  not  a  50**^  iteration,  and  where  (c)  RMP^  is  solved  otherwise. 
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Figure  22.  Master  problem  sub-optimal  solution  strategy:  Convergence.  RTS  Two 
Area  case  without  system  restoration.  This  strategy  uses  relaxed-MP  nine 
of  every  ten  iterations,  sub-optimal  integer  solutions  to  the  master  problem 
at  every  10**^  iteration  except  every  50*^  iteration  when  the  full  master 
problem  is  solved  to  optimality.  After  2,000  iterations,  we  achieve  a  65% 
gap  in  26  minutes.  In  comparison,  the  relaxed-MP  algorithm  alone  (where 
master  problems  are  solved  to  optimality  every  10**^  iteration)  achieves  a 
62%  gap  in  1  hour  and  42  minutes  in  the  same  number  of  iterations. 
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This  strategy  performs  better  than  the  original  BDA  at  early 
iterations  (gap  and  time  per  iteration  are  reduced).  However,  it  reaches  a  point 
where  the  convergence  stalls. 

E.  COMBINED  STRATEGY 

Here,  we  present  the  relaxed-MP  strategy  combined  with  the  cut  dropping 
strategy  that  keeps  the  n  most-active  cuts.  This  combined  strategy  proves  to  be 
the  best  among  all  possible  combinations.  We  show  this  combined  strategy 
applied  to  the  RTS  Two  Area  case  with  system  restoration. 

Table  4  compares  results  for  2,000  iterations  using  the  different  strategies. 
We  take  n  =  500  most  active  cuts,  and  then  combine  it  with  the  relaxed-MP 
strategy  where  the  full  MP,t  is  solved  at  every  10**^  iteration  only. 


Technique 

Original  BDA 

Relaxed  MP 

Best  500  cuts 

Combined 

Final  GAP 

16.8% 

19.1% 

48.6% 

8.3% 

CPU  time 

75h  12  m 

26m 

1h  45m 

16m 

Table  4.  Combined  strategies:  n-most  active  cuts  with  relaxed-MP.  RTS  Two  Area 
case  with  system  restoration.  This  table  shows  the  Gap  and  Time 
comparisons  for  the  most  effective  strategies.  The  results  show  over  99% 
improvement  in  time  and  a  50%  improvement  in  the  gap  achieved  when 
using  the  combined  strategy  in  lieu  of  the  original  BDA. 

In  fact,  we  can  even  obtain  better  solutions  for  this  case  by  increasing  the 
number  of  cuts  allowed  in  the  solution  strategy.  For  example,  if  we  increase  the 
number  of  most  active  cuts  to  keep  to  700,  at  the  end  of  2,000  iterations,  a  gap  of 
only  4%  remains,  but  the  time  to  complete  this  process  increases  to  20  minutes. 
We  can  also  increase  the  number  of  iterations  to  decrease  the  gap.  Figure  23 
illustrates  results  for  the  combined  strategy  with  n  =  700  when  applied  to  the  RTS 
Two  Area  case. 
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Figure  23.  Combined  strategy,  n-most  active  cuts  with  relaxed-MP.  IEEE  RTS  Two 
Area  case  without  system  restoration.  700-most  active  cuts  are  used,  and 
the  full  MP  is  solved  every  10*^^  iteration.  The  problem  reaches  a  5%  gap 
within  20  minutes  and  proves  optimality  in  45  minutes. 

The  solution  obtained  through  this  combined  strategy  represents  a 
99.12%  improvement  over  the  time  required  by  the  original  BDA.  (Comparison  is 
made  for  a  16%  gap,  which  is  the  gap  achieved  by  the  original  BDA.) 
Additionally,  this  combined  technique  allowed  us  to  prove  optimality  through 
complete  convergence. 
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VII.  CONCLUSIONS  AND  RECOMMENDATIONS  FOR  FURTHER 

RESEARCH 


This  chapter  summarizes  the  most  important  findings  of  the  thesis,  and 
proposes  areas  for  future  research. 

A.  CONCLUSIONS 

This  thesis  has  expanded  upon  the  models  and  solution  methods  of 
Salmeron  et  al.  (2003,  2004)  for  optimal  interdiction  (by  terrorists)  of  electric 
power  grids,  using  limited  offensive  resources. 

Our  first  task  was  to  validate  the  DC  power-flow  model  that  forms  the  core 
of  our  interdiction  models.  We  compare  power  flows  computed  through  the  full 
AC  power-flow  model  provided  by  the  PowerWorld  Simulator  (PowerWorld  2003) 
to  those  computed  by  our  own  DC  model,  DCCPF.  The  IEEE  Cne  Area 
Reliability  Test  Case  (IEEE  1999-1),  with  several  versions  representing  different 
interdiction  scenarios,  make  up  the  set  of  test  problems.  The  average  deviation  in 
power  flows  across  all  scenarios  is  less  than  5%,  and  all  lines  showing  deviations 
over  10%  carry  a  negligible  fraction  of  the  system’s  total  power.  These  results 
indicate  that  the  DC  power-flow  model  is  acceptable  for  interdiction  analysis. 

After  consolidating  the  mixed-integer  formulations  of  the  interdiction 
models  (with  and  without  system  restoration)  developed  by  Salmeron  et  al.,  we 
have  demonstrated  that  Benders  decomposition — this  involves  iterating  between 
a  mixed-integer  master  problem  and  a  linear-programming  subproblem — is  a 
viable  technique  for  solving  these  problems.  We  use  interdiction  decisions  as 
“complicating  variables”  and  develop  the  decomposition  methodology  through  a 
small  example  first.  We  then  extend  the  procedure  to  a  generic  power  grid,  apply 
it  to  larger  test  problems,  and  find  that  convergence  is  too  slow  for  practical  use. 
For  example,  an  instance  of  the  IEEE  Two  Area  Reliability  Test  Case  (IEEE 
1999-1),  which  has  48  buses,  requires  two-thousand  iterations  to  close  the  gap  to 
16%;  and  that  takes  75  hours  on  a  3.0  GFIz  personal  computer. 
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Not  only  is  the  number  of  iterations  large,  but  each  one  requires  us  to 
solve  a  mixed-integer  master  problem  that  becomes  more  difficult  to  solve  as  the 
iterations  proceed,  i.e.,  as  more  Benders  cuts  are  added  to  it.  For  example,  in 
the  problem  mentioned  above,  initial  master  problems  solve  in  less  than  one 
second,  but  by  iteration  1,900  they  are  taking  up  to  up  to  600  seconds  to  solve. 
(Subproblem  solution  times  remains  stable  and  short  throughout  for  this  problem 
and  all  others  tested.) 

In  order  to  improve  the  efficiency  of  the  original  decomposition  algorithm, 
we  propose  some  refinements  to  the  way  the  master  problem  is  solved.  Among 
the  techniques  investigated,  this  combination  works  best:  (a)  Solve  the  full 
master  problem  only  periodically,  say,  every  iteration,  and  otherwise  solve 
the  master-problem’s  linear-programming  relaxation,  and  (b)  drop  certain 
“unimportant”  Benders  cuts  to  limit  the  size  of  the  master  problem  (we  keep  at 
most  n  cuts,  those  that  are  estimated  in  some  way  to  be  “the  most  important.) 
Computational  times  drop  dramatically  with  this  strategy.  For  example,  a  4.7% 
gap  is  reached  in  2,000  iterations  for  the  previously  mentioned  problem,  but  this 
is  accomplished  in  a  mere  20  minutes.  Thus,  our  improved  techniques  represent 
important  steps  toward  solving  large-scale,  real-world  interdiction  problems. 


B.  RECOMMENDATIONS  FOR  FURTHER  RESEARCH 

We  identify  the  following  areas  of  research  from  which  the  existing  work 
could  greatly  benefit: 

1 .  Further  reductions  of  the  computational  burden  imposed  by  the 
master  problem.  Further  research  is  required  in  the  areas  of 
valid  inequalities  and  cut-dropping  techniques,  among  others. 

2.  Enhancements  to  speed  solutions  of  the  Benders  subproblems. 
The  current  implementation  does  not  take  advantage  of  the 
decomposable  structure  of  these  problems,  which  arises  when 
the  goal  of  the  interdictor  is  to  maximize  total  unmet  demand  for 
energy  (unmet  demand  for  power,  integrated  over  the  time)  and 
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when  repair  times  for  interdicted  components  can  vary.  This 
enhancement  is  not  important  for  small  test  problems,  but  it 
should  not  be  disregarded  for  larger,  real-word  problems. 

3.  Testing  on  real-world  problems. 

4.  Embedding  of  our  bi-level  model  into  a  tri-level  “system- 
protection  model.”  The  ultimate  goal  of  the  study  of  electric-grid 
interdiction  is  to  help  analysts  develop  plans  that  minimize  the 
potential  for  system  disruption.  Formulation  of  a  tri-level 
system-protection  model,  and  provision  of  solutions  through 
exact  and  heuristic  methods  are  thus  required.  The  system- 
defense  model  introduced  by  Israeli  (1999)  can  serve  as  the 
foundation  for  such  work. 
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APPENDICES 


APPENDIX  A:  MIP  MODEL  WITHOUT  SYSTEM  RESTORATION 
A.1  The  Interdiction  Modei  (l-DCOPF) 


This  interdiction  model,  l-DCOPF,  attempts  to  maximize  the  interdiction 
cost  (that  is,  the  sum  of  generating  costs  plus  load-shedding  penalties  provided 
by  DCOPF)  by  choosing  an  optimal  set  of  components  to  be  interdicted.  The 
following  notation  is  required  in  addition  to  that  specified  for  DCOPF: 

Additional  sets: 

C  :Interdictable  lines  (directly  or  indirectly) 

G  * :  Interdictable  generators  (directly  or  indirectly) 

Additional  model  data: 


[O  otherwise 
Jlifg.C- 
^  0  otherwise 


jlif/er 

[O  otherwise 


1  ifgG  G** 
0  otherwise 


i{g)  =  Bus  for  generator  g 


s{i)  = 


[  Substation  for  bus  i,  if  any 
[O,  otherwise 


( substation  for  bus  i  that  generator  g  is  connected  to,  if  any 
Remark:  s(i(g))  =  < 

lo,  otherwise 
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I-DCOPF; 


max  min 

SeA 


i  ge  Gj  i  c 


s.t. 


Admittance  constraints,  prior  to  linearization: 

!e/’|teLf"  seS*|/ei"  n^L’\ll<^Ll" 

which  are  re-written  below  after  linearization  of  the  50  products: 


Z  S‘“+  Z  ‘^"+  Z  <?""■)  V/ (;r,0 


iel  keif“ 


ie5'|tei“ 


//ei  l/ZGif”'’ 


Power  balance  constraints: 


/€/  i/Gir" 


//gL  |//GLf 


Z  '  -  Z  +  Z  Pt'~ +Zk = Z4  V''  ('"") 

geG,.  /|o(0=i  c  c 

Line  capacity  constraints: 


T)Line  ^  j^Line 

v/^r 

pLine  ^  pLine  ^  ^Line~^ 

V/g  L* 

^^LCap- 

nLine  ^  j^Linez-t  oBus\ 

Pi  ^Pi  (!-<>,) 

v/,//Gr,/G7f^ 

«f) 

pLine  ^  pLine ^  ^Sub's^ 

V/,5  5G;S’,/G7f' 

«) 

nLine  ^  nLine /-t  oLine\ 

Pi  ^Pi  (1-4  ) 

v/,////Gr,//G7f'' 

«) 

nLine  n^ine 

v/gr 

pLine  ^  pLine  ^  ^Line's^ 

v/g  L 

^^LCap- 

nLine  "nLine  /■%  oBus  \ 

P,  ^-Pi  (l-^^  ) 

V/,  /  /  G  r,/G  7f“ 

i<f) 

pLine  ^  pLine  ^  ^Sub  ^ 

V/,5  5G;S’,/G7f'’ 

(<) 

pLine  ^  "pLine  ^  ^Line  ^ 

V/,// //g  r,//G  zf'"' 

(<) 
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Power  generation  constraints: 


pGen  ^  pGen 
g  ~  g 

(<») 

pGen  ^  pGen /A  ^Gen\ 

g  ~  ^g  ^  ~^g  ' 

VgeG’ 

«) 

pGen  ^  ~nGen /-i  cBus\ 

^g  -^g  ^^~^i(g)' 

Vg|/(g)er 

«") 

pGen  ^  pGen  ^^_gSub  x 
g  g  ^  s(i(g))' 

VgkOXg))^o,50‘(g))G 

Upper  bounds  on  the  load  shedding: 

S.e^d,^ 

\fi,c 

/  —Load 
(^i,e 

Variable  sign  constraints: 

pGen  y  Q 

V7g  ge  <7, 

pLine  jjp^ 

V/ 

O 

Al 

\/i,c 

d:  URS 

V/ 

where 

cLine  cBus  cSub  s^Gen  ^  f  a  ill 

\d,  ,S,  ,S^  ,S^  G{0,lj| 

^  ^  X  ’  X  ’  ^^Gen  ^Gen,t  ^  X  ’  ^^Line,t  ^Line,t  X  ’  ^Bus,l  X  ’  j^Sub,t  ^Sub,t  ^  ^ 

geG*  leJ^  ief  seS* 

A.2  The  Dual  Interdiction  Model,  DI-DCOPF 

Taking  the  dual  of  the  inner  model  presented  in  section  A.1  yields  the  dual 
interdiction  model  shown  below. 
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DI-DCOPF: 


'A^sr  «■-<•)+  X  sf'  ) 


hi\i€l  ,/Glf 


+  Z  <?■*  «-^/*)+  Z  «-<*) 


/,//|//Gi  JleLY 


+Z<"  <Z4)+Z^’X‘  +  Z/i°(i-‘?")<+  Z 


g|‘(g)e^ 


gk{'(g))e5 


/,?|/g/  ,/GZ-f' 


hs\seS  JeLf 


+  X  /;“'a-<5;r)(<,f-<)+Z4  < 


l,ll\lleL  Jleiy 


Power-generation  dual  constraints: 

<) + (1  -  )<“  +  ^  K 

Line-flow  dual  constraints: 

Tuf  +  Tuf  -I-  (1  -  /l;^  -I-  7rf°  )  -I-  A.I'TU^  +  A.l'Tuf  -I-  ^  +  Tuff  j 

/,ijfG/*,/GZ,. 

,  /  ^LS'  ,  ^LS'^  \  .  ( ^LL  ,  ^Lt  \  ^Bal  ,  ^Bal  n 

+  L  +^i,s  )+  2.  [^m )-^oii)+^dii)=^ 

;  ,s  IsE  S*  ,Z6  i,  Z  ,z/|zze  i*  ,/Z6  if 


Power-shedding  dual  constraints: 


—Bal  ,  —Load  ^  r 

^fic 


V  i,c 


Phase-angle  dual  constraints: 

-  Z  )  +  Z  ) = 0 


Signs  on  dual  variables: 
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unrestricted 


,—  -A.  ,,m^Lf\  -—Z 

n 

,77 

.71 

<0 

-I7,  L* 

71  ,71^,71  , 

n 

-pY 

.71 

>0 

„G  —Go  — GS 

77  ,7t  ,77  . 

1  71 

<0 

^Load  <  Q 


Bounds  on  the  dual  variables: 

TT  \7r  ,7r  ,K  -.-K  = -Max.  shedding  cost -Max.  Generating  cost 
7t^  ,7t^  ,7t^^  ,7t^^  :-7i^  =-Max.  shedding  cost -Max.  Generating  cost 

7t^\7t^  ,71^^  ,7t^^  ,7t^^  :  .^^  =  Max.  shedding  cost  +  Max.  Generating  cost 
71^^  ■.  -Bfi 

71^*  : 

where  we  take  9  =  \  radian. 

A.3  Linearization  of  DI-DCOPF 

Linearizing  the  57t  products  in  DI-DCOPF  yields: 


LDI-DCOPF: 


max^M,]  -v/")+  ^  X  -1^/^)+  Z 


,/Glf 


Us\s^S* 


/,//|//gL  ,//Gi;‘ 


s\s(i(s))£S 


+ z  (^1°  -  ) + z  )  -  (vf  -  vf )] 

l€  I^*  /e  Z,’^ 

I  X  -j-\Line  /  —LB  —LB^  \  /  LB  LB^  \  .  X  ’  j-iLine  /  —LS  —LS^  \  /  ZS  LS^  \ 

+  L  Pi  -^i,i  )-(yi,i  )+  L  Pi  -^i,s  -vi,s  ) 

l,i\ief  ,leLf“  Z,s|i65*,/6i“ 

+  X  <*'+ZZ4  4“" 

.  ..  1—  — I  • 


/,//|//gZ  ,//GZf 


liR)  ^  ^  s  R  R  i{^Y  p:  s(i(2)Y*'g:  —  "'2 


Vg  (Z/) 
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Tif  +  71  f  +  (1  -  +  7t^ )  +  Xi7i\  +  X’^Ttf  +  ^  [Trf^  +  Ttff  j 

l,i\^ief  ,1&L^ 

,  /  ^LS'  ,  ^LS''  \  ,  /  ^LL'  ,  \  ,  ^Bal  n 

+  L  1^/,.+^/,.)+  2.  +^/.«  j-^o(/)+^^(/)  =  o 

Z,s  .seS*  ,/64  /,;/  ZZei*  Jlelf" 

V/ 

(P,A 

^Bal  ,  —Load  ^  r 

+^i,c  ^fic 

V  i,c 

(s,,A 

-  Z  )  +  Z  [^f 

l\i=o(l)  l\i=d{l) 

+  77f^  =  0 

V/ 

m 

Linearizing  constraints  are: 

^  —A  cLine 

<  Tt^ 

V/g  L* 

(K) 

o 

VI 

+ 

1 

+ 

> 

V/g  L* 

(rC) 

A^  -v.  ;=^ /I  cLine\ 

V,  -77,  >-7r,{\-d,  ) 

v/g  L* 

(<) 

BA^  ^  —A  cBus 

v/,/|/Gr,/G7f“ 

(O 

o 

VI 

+ 

1 

+ 

-2 

v/,/|/Gr,/G7f“ 

(^;) 

BA'^  —A'^  -s.  ;=A/i  cBus\ 

^14  -^l  ) 

v/,/|/Gr,/G7f"' 

(^1) 

SA'^  ^  —A  cSub 

VsJlseS^JeL^;'’ 

(O 

O 

VI 

+ 

1 

+ 

■^'‘ 

VsJlseS^JeL^;'’ 

(Ch) 

SA'^  —A'^  ;=A  /-I  cSub  \ 

^l,s  -^l  (1-^^,  ) 

'^s,l\sGS\leL^;'’ 

(^>3) 

LA^  ^  ■=  A  QLine 

V/,//  ^  \ 

yi,ll\lleL\lleLA' 

<. 

1=:  ^ 

+ 

I 

+ 

lA 

o 

V/,//|//Gr,//G7f'' 

(^ih) 

LA^  —A^  -v.  ;=A  /I  Saline  \ 

\ii -^i  ) 

V/,//|//Gr,//G7f“'' 

(7(1,11%) 

LA  ^  —L  s^Line 

V,  <  77,  5, 

V/g  L* 

(K) 

o 

VI 

^  ■,. 

1 

t]  ^ 

V/g  L* 

irt) 

vf  -77 f  >-^l^{\-dA"A 

V/g  L* 

(iA:) 
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V 


LB^ 

14 


LB^ 

V;,- 


LS^ 

LS^ 

LS^ 
V,  . 


LL^ 

141 

Lt 

; 

141 

Lt 

; 

141 


V, 


'^i,i 


BA^ 


SA- 

SA- 

V,.. 


V 


SA- 

l,s 


LA- 

Vl,ll 

LA- 

\ll 

LA- 

V,J, 


^  =i  cBus 
<711  Oi 

v/,/|/Gr,/Gif"" 

(rSO 

-<r<o 

V/,/| /£/’,/£  If"' 

(^1) 

-K’^f 

V/,/|/Gr,/Glf"^ 

(^1) 

< 

'^s4\sgs\Igl^;'’ 

0 

VI 

+ 

<>2 

1 

V5,/|5G5*,/Glf' 

—LS^  >>  CK  S^Sub  \ 

V5,/|5G5*,/Glf'’ 

V/,//|//e  r,//e  If^ 

i7(i,u\) 

0 

VI 

1 

V/,//|//e  L\IIg  If^ 

(r^LuO 

-<  >-^;(i-C') 

V/,//|//e  L\lle  If^ 

(r^uJ 

V.  :=A  cLine 

>  -TT,  Si 

iK'> 

1 

IV 

0 

V/g  L* 

(K) 

-Trf 

V/g  L* 

(K^ 

> 

V/,/|/Gr,/Glf"" 

(C) 

1 

IV 

0 

Vz,/|/Gr,/Glf"^ 

v/,/|/Gr,/Gif"' 

V5,/|5G5*,/Glf'’ 

1 

IV 

0 

V5,/|5G5*,/Glf'’ 

V5,/|5G5*,/Glf'’ 

:=^  cLine 

>  -TT,  S„ 

v/,//|//Gr,//Gif"'' 

ifihii),) 

1 

IV 

0 

v/,//|//Gr,//Gif"'' 

-Trf  <<(l-4f“) 

77 


vf  > 

V/g  L* 

(^) 

v/"  -;r/'  >0 

V/g  L* 

(r^;) 

1 

VI 

1 

^  ~~ 

> 

v/g  L* 

(rl;) 

LB  -s.  :=L  cBus 

v/,/|/Gr,/G7f“ 

(C) 

LB  ^—LB  'v.  A 

-^i,i  ^0 

v/,/|/Gr,/G7f“ 

(7(1, Oi) 

LB~  ^LB~  ^  —L  /I  cBus  \ 

-^14  ) 

V/,/|/G/,/G7f“ 

(7(1, ih) 

1 

Al 

U  <>5 
•>J  ^‘' 
> 

V5,/|5G5*,/G7f'’ 

(?^').) 

7,5*  ^BS  A 

-^l,s  ^0 

VsJlseS^JeL^;” 

LS~  ^BS~  ^  ;5F^ /I 

V,,s  -^l,s  ) 

VsJlseS^JeL^;” 

(7^1%,) 

>  -Ttfdr 

V/J/IZ/Gr^/Glf'"' 

(7(1,11),) 

V/,//|//g  L*,//g 

(T^uih) 

\n  -Kfl  ^^"(1-4) 

V/,//|//g  L*,//g 

(T^uoJ 

G  V.  cG 

V  >  —K  O 

g  g  g 

VgGG’ 

(7D 

<-<^0 

VgGG* 

(7D 

VgGG’ 

(^J 

GB  V.  —G  cBus 

V  >  —7t  O  t  . 

g  g  ‘(g) 

Vg|/(g)er 

(7D 

vf -;rf  >0 

Vg|/(g)er 

(7::) 

-.,G5  ^GB  ^  =0/1  J^Bus\ 

Vg|/(g)G  I* 

(7D 

Vg|^(/(g))G 

(7D 

vf-<^>0 

Vg|^(/(g))G 

(7^:) 

^,GS  ^GS  ^  ^0/1  cSub  \ 

^g  -^g  ^^g(^-^sii(,))) 

Vg|5(/(g))G 

(7::) 

7cf  <nf 

V/ 

(vf) 

Ttf  <  TTi 

V/ 

(vf) 
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^LB*  ^  —L 
^l,i 

V/,/|/G/*,/G7f" 

(nu) 

^LS*  ^  —L 
^l,s  - 

V/,5|5G  S\l^  Lf’’ 

(vtf) 

Kn 

V/,//|//G7*,//G7f“'' 

ivf!u) 

IV 

1 

V/g  L* 

inf) 

IV 

1 

-Z^l 

V/g  L* 

ivf) 

^LB~  ^  —L 

Kii  >-77i 

V/,/|/G/*,/G7f" 

ivff) 

IV 

1 

-r^i 

t~- 

V/,5|5G  S*,lG  if* 

ivff) 

it=^ 

1 

Al 

^  =5 
“-1 

V/,//|// Gif //Gif'" 

ivff) 

VgGG’ 

i€) 

Vg|/(g)er 

ivT) 

Vg|5(/(g))G5* 

ivf) 

71'^“^  unrestricted 

^LS  ^LL  A 

7U  ,K  ^  ,K  ,K  ,7Z  ,7t  >0 


>0 


jmm^-  >i«v^n  /~\ 

7t  ,7t  "  ,71  ,  Tt  ,K  ,K  <0 

,  v"" ,  <  0 

7t^  ,7t''\7t^\  <0 

^Load  <  Q 


A.4  The  Master  Problem 


The  master  problem  for  Benders  decomposition  derives  from  the  dual  of 
LDI-DCOPF.  At  iteration  k,  the  master  problem  has  the  form: 

(MP^.):  max  z 

<5‘€A,zg]R 
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s.t. 


leL 

;,!|fe/*,/eif“' 

/,5|.VG5",/GLf^ 

,  X  ^  saline /=  A  .±A^  ,  =z4-Xz4‘^  =A.±A~  :=A.±A~  \ 

+  2^  4  (^1  7(!j!)„k'+^i  7iij!)„k'~^i  7iiji)„k'~^i  7(iM)„k') 

!,ll\IleL\lleLf‘^ 

+T.^r(^"rl'ik'+Kr^ik'-Kri:,k'-KK,k') 

tei* 

+  ^  ) 

^  ‘■^f'  f(l,i\,k' ^ '‘‘l  f  (l,i\,k'  '^l/(l,i\,k'  '^l  /  {l,i\,k') 

l,i\iel\leLf“ 

^  ^  V‘'l  /{l,s\,k'^  ^‘'l  /(l,s),,k'  ^‘'l  /(l,s\,k'  •'^l  /(l,s),,k'} 

l,s\seS\leLf‘’ 

I  cLine  /  yLL  ,  yLL  _  yLL  _  =L  \ 

IMllel'jleLf" 

glgeC  slgeG 


+  z.c,(-^c-^;C')+  Y.Kf:r 

S\i(sYl  S\i(gYi 

+  Z  ST\-Ky..-Kfiy  Z 

g|i(!'(g))e5  g|s(i(g))eS 

l\l  6i*  lAiefjelf"  Z,s|i6S*,tei“ 

!,ll\lleL' JleLf"  1\1  el'  !,i\iel’ JeLf“ 

■*■  Z  ('y{i,sh,k' ~7(i,sh,k')'^  Z  ('y{i,ii%,k' ~  7{i,ii'),,k') 

l,s\seS\te  i“  /,Z/|Zte  L*  ,ZZe  if" 

+zr(^A^i -^^'))+(^z"(75. -<i'))i+  z  r(^/^(^(w' -Cu')) 

I  *-  ;,i|i£/,/sif“ '- 

+  Z  r(^"(Cu'-C),r))l+  Z  \[^UV(ui),k’-vZi)r)) 


l,s\ssS  ,/e£7 


/,//|//eZ  ,ll^L, 


+Z(-»f*'>+ Z  (-»f’z">+  Z  (-»f<r)+Z*./f:.'+ZAA 


(?,c),^' 


geG 


<{g)e^ 


s('(g))eS 


V£'  =  l,2,..l 
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APPENDIX  B:  MIP  MODEL  WITH  SYSTEM  RESTORATION 


B.1  Constructing  Time  Periods,  and  Additionai  Notation 

The  notation  below  and  time  period  construct  algorithm  are  extracted  from 
Salmeron  et  al.  (2004). 

Required  notation: 

T  =  set  of  periods,  for  T 

^  =  Z*  u  G*  u  5*  u  5** ,  set  of  all  (directly)  interdictable  elements 
Dur(e)  =  Duration  (hours)  of  outage  for  element  ee  ^ ,  if  attacked 
=  Duration  (hours)  of  time  period  t,  for  re  7 

{1,  if  component  e  remains  unrepaired  in  time  period  t  after  being  attacked 

, 

0,  if  component  e  is  repaired  before  time  period  t  after  being  attacked 
for  teT ,  ee  ^ . 

Remark:  In  the  above  notation  pf;; ,  and  pff  denote 

when  e=l  is  a  line,  or  e=i  is  a  bus,  or  e=g  is  a  generator,  or  e=s  is  a  substation, 
respectively. 


The  following  algorithm  constructs  the  set  of  time  periods,  T,  based  on  the 
different  outage  durations  for  all  interdictable  elements.  In  the  course  of  the 
algorithm,  and  are  also  constructed: 
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Algorithm  “Construct  Time  Periods”: 

Initialization:  ^  T  <—  0,  t  <—  0,  <—  0; 

While  ^  ^0: 
t  < —  r  + 1 

[O,  otherwise 
m,  <—  min|£)Mr(e)|eG 

^  ADur{e)  = 

Z), 

End  While 


Additional  notation  for  interdiction  model: 

I^*  =  Set  of  lines  /  that  ean  be  direetly  or  indireetly  interdieted  in  period  t. 
Note  ; 

/  G  e7  if  either: 

=  or 

y5®"^  =  lforsome/|/GEf'“,  or 
fixT  -  i  for  some  5  I  /  G  Ef  *,  or 
y57i“  =  lfor  some // 1 //Gif'- 
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|l  if/G  L* 

[O  otherwise 

|i  if/Gi: 

[O  otherwise 

|l  ifgeG’ 

0  otherwise 


If 

[O  otherwise 

a;  =  \^  if^G/ 

0  otherwise 


|l  if^G^’ 

[O  otherwise 


B.2  The  Interdiction  Modei,  l-DCOPF-R 


l-DCOPF-R: 


max  min 

s  (p<^‘\p^‘-‘,s,A)  ,=i  [  s  ,  c 


S.t. 


iel  lleLT 


seS  \leLv‘ 


ZoLine  cLine  \ 

Pm  Pi  ) 


V/,Vt 


«i) 


lie  nils  lY'"" 


P^r *  I  Pf:-S‘-+  X  ^“<5,“  + 


iel  \leLr^ 


seSlleLf^ 


ZoLine  cLine\ 

Pm  4  ) 


V/,Vt 


lleL  \lleLr 


i^r-  z  p!T‘+  z  ^;"'+z^«,=z<. 

geGi  l\o{l)=i  l\d{l)=i  c  c 


«,) 

«f) 
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jyLine  ^  -pLine 
^t,l  — 

(^5) 

pLine  ^  pLine  ^  ^ Line's^ 

'^tje  L\  I5^f^=\ 

«/) 

riLine  ^  T^Line/i  cBus\ 

P,j  ^Pi  (1-^/  ) 

pLine  ^  pLine  ^  ^Sub's^ 

Vt,/,5  5G5*,/G4,y5“=l 

pLine  ^  pLine  ^  ^ Line's^ 

Vt,/,////Gr,//Glf'',y5,^"=l 

K]i) 

pLine  "pLine 

«?) 

pLine  ^  pLine  ^  ^  Line's^ 

«/) 

j~%Line  riLine/-t  cBus\ 

P,j  ^-Pi  ) 

vt,i,i\i&i\  =  i 

nLine  riLine/-t  cSub\ 

Ptj  ^-Pi  (1-^.  ) 

«0 

pLine  >  _pLine^^  _  ^Line^ 

Vt, /,//// G  r,//G  If 

«S/) 

pGen  ^  pGen 
^Lg  -  ^g 

VI  g: 

(<D 

pGen  ^  pGen  /i  ^Gen  \ 

^Lg  ~  ^g  ^g  ' 

Vlg|gGG*,AS”=l 

«) 

pGen  ^  "pGen  /i  cBus  \ 

^Lg  -^g  ^^~^i{g)) 

vig|/(g)Gr,A.ff=i 

«) 

pGen  ^  pGen  /i  cSub  \ 

^Lg  -  ^g  ^^~^s(i(g))) 

Vlg|.0'(g))G5*,Aff,„  =  l 

«') 

S,,,e  ^  d,c 

\fi,c 

/  —Load  \ 

) 

pGen  y  Q 

'^t,g 

pLine  jjpg 

Vt,l 

SiMe  ^  0 

\/t,i,c 

^i,i  URS 

Vt,i 

Bounds  on  : 

7r‘^\7r‘^ =  -Max.  shedding  cost -Max.  Generating  cost 

,7t^^  ,7t^^  -.-Ti^  =-Max.  shedding  cost -Max.  Generating  cost 
,7r^  ,n^'^  ,7t^^  ,7t^^  :  .^^  =  Max.  shedding  cost  +  Max.  Generating  cost 
■.  -7i^  -Bfi 

where  we  take  9  =  \  radian. 
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B.3  The  Dual  Interdiction  Model,  DI-DCOPF-R 


Taking  the  dual  of  the  above  primal  problem  yields: 
DI-DCOPF-R: 


max  max  y  ivi,  y  < 

5 


«-<)+  z  «-<) 


iel  heir' 


+ 


+  Z  «-<,*)+  Z  «-<)■ 

seS’|/eif'’  lleL‘\leLf;" 

Z(Z4)  <'+ZZ^>5+Z  Z 

i,t  c  t  gtG”  geG’t\p°l"=\ 

+  z.  z  ^;(i-cx+  z  .  z 

s\i{g)el  «IA?('j;)=l  gk(‘(g))e-5  t\Pf%(g)) 

+ZZ^‘(^5-^5>+Z  Z  >(<;-<;> 

i  tei**  /ei’ 

+  Z  Z  Z  Z 

le  /*  |/e  if “  =1  seS''\le  if*  r|/3“  =1 

z  z  ^'(i-4r)(^“.-<;,)+Z‘'.. 

lie i* |/e 4"  i| A;,"  =1 


-I- 


s.t. 


+(1-A4^S  +4< s V  <,g  (/;“) 


+ (1  - 1,:;  + 4 ) + +  z  + < ) 

i6/'|Zeif“,/?,4’=l 

+  z  (</XS)+  2^  (<//+<//)-<;/)+<i)=o 


seS'|tei“„5“=l 


ZZ6r|ZZe4",/?,^r=l 


V4/,c(S_) 
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-  Z  [Ki  +  )  +  Z  ) = 0 


^t,i  (0.i) 


;r^“^  unrestricted 
^  ,7t  ,7t  ,  7t  ,7t  ,7t  >U 

7t^  ,7t^\7t^  ,  ,7t^^  ,7t^^  <0 

_Gr,  ^GB  ^GS  ^  r\ 

7t  ^  ,7t  ,7t  ,  7t  <  U 

<  0 

B.4  Linearization  of  DI-DCOPF-R 


Linearizing  the  products  Stt  yields: 
LDI-DCOPF-R: 

"s  L  pLine  /■A'  \  ■  X  ^  pBus  /  BA'  BA^  \ 

2^  Pui  iytj  -\i  ) 


max 

S.Tt.V 


iel  lleLT 


lleL'lleLfr 


seS  lleLv' 


+Z(Z4)-^,“+II^>5+Z  Z 


id  c 

t  geo" 

geG’  I 

•\P?T=^ 

+  z. 

z 

Z 

g|<(g)er 

g\s(,(g))eS 

+ZZ 

S)+Z  z 

p,"  («; 

„L* 

^t,l 

t  l€L** 

tei*  t\fi,f‘  = 

1 

+  Z 

Z 

'„LB-  „LB*  \ 

/  LB~ 

LB*  \\ 

\l,i  )J 

fe/*/eif“ 

'  t\pf“=\ 

+  Z 

I.  PH 

/„i5‘  \ 

/  LS~ 

LS*  H 

+  Z 

Z  PH 

( 

lleL'\leL‘{r 

J.GS  GS) 

'^t,g  ''t,g ) 


Load 


die 


With  balance  constraints: 

<i) + (1  -  Ki  )<i + 4"< + 4.)  ^  D{t)h 


(C) 
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5)  +  l‘/3"”< 

,  0  Z,  oLine  —L' 

+  ^iP,j 

■  +  Z 

1  \ 

/ 

+ 

■  Z  ( 

—LL  ,  _ii+  \ 
^t,l,ll  +^tj.ll  j 

—Bal  ,  ^Bal  A 

V  t,l 

iPu) 

/tei* //eif" 

^Bal 

V  t,i. 

^  (S,,c) 

1 

M 

^/(<  +<")+Z^/( 

<  +^^)  =  o 

Vt,i 

(^u) 

l\o(l) 

/i^/(/) 

Linearizing  constraints: 


^  ^  Saline 

\/ 

'^t,i\i&i:,t&T,pA =\ 

(r^,iA 

<;-<<o 

\/t,l\lGL\t&T,/3A‘'  =  l 

(r^jA 

\/t,l\l&  L\t&T,/3A"  =l 

(rAA 

yt,i,l\i^  f  ,1^  LA  ,t^T,PA=^ 

(C).) 

BA*  „A*  ^  r. 

\l,i  -^ui  ^0 

(AS, A 

BA^  —A^'  -s.  7=A /■■%  cBus  \ 

\i,i -^tj  ) 

'\ft,i,l\i^l\l^LA,t&T,pA  =  '^ 

(AS, A 

SA'^  ^  —A  cSub 

\l,S  ^  ^tAs 

'^t,s,l\s^S\l^LA,t^T,PA  =  ^ 

(AS,  A 

'^t,s,l\s&S\l^LA,t^T,l3A  =  l 

(AS,  A 

SA'^  —A'^  v.  —A  /I  S^Sub\ 

\l,s-^t,l  ) 

\ft,s,l\s&S\l&LA,t^T,PA  =  ^ 

(AS,  A 

LA^  ^  —A  oLine 

\ui  ^  ^tAii 

vr,/,//|//G  r,//G  LA,t& T,pA'=\ 

(AS,  A 

LA*  ^A*  ^  n 

\lM~Al  -0 

vr,/,//|//G  r,//G  LA,t& T,i5A"=\ 

(AS,iA 

LA^  —A'^  v.  —A  /I  cLine\ 

vr,/,//|//G  r,//G  LA,t& T,pA'=\ 

(AS,iA 

lA  ^  —L  oLine 

\i  ^  ^tAi 

\ft,l&L\t&T,PA=l 

(A,  A 

\l  -^ui  ^0 

(AA 

(AA 
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LB^  ^  =Z,  S^Bus 

yt,i,i\i& r,ie  LA,t& r,y5,7=i 

^rAlA 

LB*  „LB*  ^  c\ 

^tu  -^,j,i  ^0 

yt,i,i\i€i\i€LA,t^T,j3A='^ 

ATlA 

LB^  ^LB^  v.  —L  /I  s^Bus  \ 

'^t,i,i\i&r,i&LA,t&T,/3A=^ 

(^1)3 ) 

LS'^  ^  =1  cSub 

-  ^tAs 

Vt,s,l\s€S\l€LA,t^T,Af='^ 

^t,l,s  -^LUs  ^0 

^t,s,l\s€S\leLA,t^T,Af='^ 

(AS,sA 

LS^  —LS^  >s.  -=1  /I  ^Sub\ 

'^t,s,l\sGS\lGLA,t^T,/3A='^ 

(AS, -A 

LA  ^  ;=L  cLine 

\lM  ^  ^,All 

yt,i,ii\ii& r,iiG  Lf  ,^G  t,at=^ 

(AhA 

lA  _z,z,^  ^  (\ 

^t,i,ii  ^t,i,ii  —  0 

r,iiG  Lf  ,^G  t,at=^ 

At  J, 11)2  ^ 

LA  _ LA  \  —L  /i  s^Line\ 

V^,/,//|//g  r,llG  Lf  ,^G  T,AT=^ 

At, 1,11)3^ 

A~  —A  cLine 

\i  ^-^tAi 

yt,i\i^i:,t^T,pA=\ 

Aa"> 

<-<>0 

(AA 

A  a) 

Au  ^ 

\ft,i,l\isl\l^LA,t^T,PA=\ 

AtlA 

BA  /rr-^  A 

-0 

'^t,i,l\isl\l&LA,t&T,j3A=\ 

ATa"> 

BA~  ^A~  ^  :=A  /I  S^Bus  \ 

\i,i -^,j  ) 

yt,i,i\iGi\i&LA,t&T,j3A=^ 

AT, A 

AL^-KAA 

'^t,s,l\sGS\lGAA,tGT,/3A='^ 

AT)) 

AL-<  ^0 

^t,s,l\s€S\l€LA,t^T,Af='^ 

ALa"> 

SA~  —A  ^  =A  /-I  cSub\ 

\l,s-^.J  ) 

^t,s,l\s€S\l€LA,t^T,Af='^ 

AaA 

LA~  ■=A  saline 

v,jj,  ^  -^,Aii 

v^,/,//|//g  r,//G  LA,t& t,at=^ 

AT  A 

LA  ^—A  'v.  A 

v,j,ii-^tj  ^0 

V^, l,ll\llGL\llGLA,t^  T,  at  =  1 

AT  A 

LA~  ^A~  ^  —A  /i  cLine\ 

v^,  /,  // 1  //  G  r ,  //  G  ^  G  T,  AT = 1 

ATa) 

v^;  >  -AaA^ 

V^,/|/Gr,tGT,y5,^”‘=  =  l 

A, 1)7 

Ai-Aj  ^0 

v^,/|/Gr,tGT4y"^=i 

At,i)i^ 

v^,/|/Gr,tGT,y5,y“=i 

Ala) 
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15  V.  —L  cBus 

V^,/,/|/Gr,/Glf%tGT,y5,7*  =  l 

At,i,A 

yt,i,i\i&  i\i&  A\t&T,Ar='^ 

AmA 

LB~  ^LB~  ^  =Z,  /-I  S^Bus  \ 

yt,i,i\i&  i\i&  A\t&T,Ar='^ 

AlA 

LS~  \  —L  cSub 
\l.s  ^  -^,As 

yt,s,i\s&s\i&LA,t&T,pA 

Am  A 

z,s  n 

\l,s  ^0 

yt,s,l\s&S\l&LA,t^T,PA  =\ 

AAA 

'^t,s,l\s&S\l&LA,t&T,pA  =\ 

AAA 

^  <AA^ 

v^,/,//|//g  r,//G  LA,t& T,AA=^ 

AAA 

Am  -AA  ^  0 

yt,i,ii\ii& L\u& LA,t& t,aT  =  ^ 

AAA 

Aui-AA^Ka-A) 

v^,/,//|//g  r,//G  LA,t& t,at=^ 

AAA 

<A<As 

yt,g\g^G\t^T,pA=\ 

A  A 

\ft,g\g^G\t^T,PA=^ 

A  A 

yt,g\g^G\t^T,/5A=\ 

A  A 

vS^-CC) 

^t,g\i{g)^l\t^T,PA,)=^ 

AlA 

^t,g\i{g)^l\t^T,PA,)=^ 

AlA 

GB  ^GB  ^  :=G  /-I  cBus  \ 

yt,g\i{g)^i\t^T,p:-^^=\ 

AlA 

G5  :rrG  s^Sub 

\g  ^ -AAs(Ks)) 

yt,g\s{i{g))^s\t^T,pfA,))=^ 

AA) 

yt,g\s{i{g))^s\t^T,pAM,),=^ 

AlA 

^t,g\s{i{g))^S\t^T,PfA,))=^ 

AlA 

^-4+  ^  =A 
^t,l  -^tj 

'^t,l\l€  L\tG  T 

A ) 

</ 

\ft,l\lGL\teT,Ar  =  '^ 

Aa 

—LB*  ^  =L 
^t,u  ^^,,1 

\ft,l,i\i^  I*J& LA,t& TAfr 

aaa 

^—LS^  ^  ^f:f:L 

^,M.  ^  Al 

\Jt,l,s\s^A  ,1^lA  ,t^TAA  =  ^ 

AAA 

—LL*  ^  =L 
^t,ljl  - 

\lt,l,ll  1  //g  L*,IIg  LA,t& TA'A  =  1 

Amu) 
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^t,i  -  ^t,i 

yt,l\l^  L\t& T 

infd ) 

'7s 

1 

Al 

'\ft,l\l^L\t^T,l37;^  =  \ 

inh) 

^—LB 

yt,l,i\iGr,lGL^‘‘\teT,ff;^  =l 

(O 

1 

Al 

yt,l,s\sGS\l€L7‘\t€T,j37f  =\ 

(O 

^LLT  ^  —L 
^t,iji  -  ^t,i 

\ft,  /,  // 1  //  G  r ,  //  G  ^  G  T,  =  1 

v^,g|gGG*,rGr,AS"=i 

(O 

v^,g|/(g)Gr,^GT,AS;)=i 

(O 

v^,g|50-(g))G5’,rGr,A;;-^,  =  i 

(O 

n^‘''  unrestricted 

V^\v^Cap^, 

,  ;r“^ ,  <  0 

^A-  ^  ^LCap-  ^  ^LB-  ^  ^1$-  ^  ^LU  ^  ^G  ^  ^GB  ^  ^G5  <  g 

7t^ ,7c'^\7t^\  7t‘^^  <0 

^Load  <  Q 

B.5  The  Master  Problem 

The  master  problem  of  the  Benders  decomposition  results  from  the  dual  of 
LDI-DCOPF-R.  At  iteration  k,  the  master  problem  has  the  form: 

(MP^.):  max  z 

J€A,zg]R 
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s.t. 


IL  +Kr(t,i),,k' -Kr^,i\,k’ -Kru,i)^,k-) 


t\P“r=\  leL 


t\/3f"-^=ll,i\teI*JeLf"^ 


t\fi^'‘^'’=U,s\seS*JeL^/^ 


I  X  ^  X  ^  Saline /■■=  A  ,  =A  :=A  ^LA  :=A  ^JA  \ 

■*"  ^  ^  ^ii  i^tjy{t,iji\,k'~^^tj7it,iji)^,k'  ^tj7(tj,ii)i,k'  ^tj7^i,u)^,k') 

t\/}kir=l  /,//|/tei*  JleLf‘‘’' 

t\Pl-f‘=\  Zei* 

+  ^  ^  /  J^L  yLB  _|_  —L  yLB  _  —L  yLB  _  —L  yLB  \ 

Z|/S*“=l  l,i\iel‘jeLf“-' 

+  ^  ^  ^Sub  /  —L  _|_  =i  yLS  _  j=L  _  =L  \ 

Zj  Z-l  y‘'t,l/ {t,l,s\,k' ^  ^‘'t,lf(t,l,s)i,k'  (t,l,s\,k'  ‘'^t,l/(t,I,s),,k'k 

r|^“=l/,i|seS' ,/€/,“ 

I  ^  ^  cLme  /  =L  yLL  ,  =L  yLL  _  yLL  _  yLL  \ 

t\Pl‘)l‘=\  Z,ZZ|ZZeL*,ZZeLf“'' 

+  E  E  E  E  CrJ 

zlA'5"=lglzreG'  z|/3,‘J"=lg 

+  z  z  c;(-<C).,r-<C)3.^')+  z  z 

‘\PuA^  g|Z{g)e/*  Z?l'(g)eZ* 

+  E  E  ^“(-KjZ>,r-KX%,r'>+  E  E 


G 

h‘,gh,k' 


^G  ,,GB 

t.gh.k' 


^IA,j;?,(s))=lgkO'(g))eS 


'IA.:?,(f;))=lgk(Z(g))eS 


+  Z  T  ^MGih,k'-r^jh,k''>+ T  Z  ^z’^(^l)3.r-^Jz)3.r) 


+  Z  Z  ^z^(^l.)3.^'-:^l)3,r)+  Z  Z 


t\^;f=lI,s\seS  ,leL, 


+  Z  Z  Z  Ki(r^i,h,k'-r^u)„k^ 

t\P^j“=l  l\l  ei*  z|AT'=1  /,!>'£/' ,Zeif“ 

■*■  Z  Z  ^t,l(7{t,l,s)^,k'  ~  7(tJ,s),,k')'^  Z  Z  ^t,l(T{t,l,ll),,k'  ~  7{t,l,ll),,k') 

z|/3“=l  Z,«|ieS*,ZeL“  l\PlT=^  UU\U^P PA" 

+ZZ^z’^(^(MU'-^(lo,^’)+  Z  Z^z!z(7(m),.’-<zu')+  Z  Z  Ki^nZ,\k'-nZ,\k^ 

t<^T  leP  t\plf‘=\l(LL  t\pf“’=\lAiel'' 

+  Z  Z  Z  Z  KMw-Vi!:i,,i),k') 

/,s|ie5’,ZeL“  t\Pl'T=^  l,tt\tteL‘ JlGlf" 

+  Z  Z(-<C)A')+  Z  Z  (-<Cu')+  Z  Z  (-<Cu') 

<;=lgeG*  r|^^;;',,=li(g)e/* 

+Zl0JrX.>r+i:i,0,f.A,uu' 

teT  g  teT  i,c 


'IA“,(3;))=l^(<{i?))e5' 


yk'  =  \,2,...,k 
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APPENDIX  C:  LINEARIZATION  OF  CROSS-PRODUCTS 


This  is  an  explicit  list  of  constraints  for  the  three-bus  example  in  Chapter 
IV.  Y,ri  represent  the  dual  variables  for  the  linearizing  constraints. 


’12  ^\1  ^\1  —  ^ 

irtlx) 

0 

VI 

1 

V  2 

iiii) 

12  ^^^12  ^12  — 

1 

u> 

lA 

0 

(A\) 

0 

VI 

1 

V  2 

ir^i) 

23  ~  -^23  ^23  -  0 

(?4i) 

1 

lA 

0 

{7^,2) 

23  ^23  "^23  ^23  —  "^23 

(2^33) 

''2+^12^12^0 

ir^n) 

''2  -<2  ^0 

'12  ~^\i  +'^12  ^12  — 

(7^22) 

'l3  +■^13  ^13  —  0 

i7^2x) 

1 

IV 

0 

(7^2) 

'l3  '^13  +  ■^13  ^13  —  ■^13 

{7^22) 

'23  +  ■^23  ^23  —  0 

(2^31) 

'23- <3  ^0 

(7^22) 

'23  ~  '^23  +  '^23  ^23  —  ■^23 

(7^22) 

LCap^  =LCap^  ^ 

12  ^12  —  ^ 

(7^2?^^) 

LCap"^  —LCap^  ^  ^ 

12  '^12  —  ^ 

(7^22^^) 

LCap^  —LCap^  —LCap^^L  ^  =LCap^ 

12  '^\2  '^\2  ^12  —  ^^^12 

(7^2^^^) 
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LCap  _^LCap  <  A 
M3  ^13  —  ^ 


LCap  _  —LCap  ^  ^ 
M3  ^ 


LCap^ 

M3 


.LCap^ 


13 


-TT, 


LCap'' 


>  -W^Cap* 

c’la  _  /t.j3 


(^3 


Cap^ 

131 


^  yLCap^ 
^  Y^Cap^ 


v^'  -3tS‘^  Si,  <0 


v,7-'  -3-““'’  <0 


LCap 

■^23 


LCap''' 


- _  TT-'^C’iip  _  7=LCap  cL  ^  _=LCap 

^23  ^^^23  ^^^23  ^23  “  ^^^23 


^yLCap^ 

ifii 


Cap' 

232 


^  yLCap^ 


LCap-  =LCap-  gL  >  q 

M2  ^^^^12  ^12—^ 


itir) 


LCap 

M2 

Leap¬ 

'll! 


-K, 


LCap 


12 


>0 


„LCap  ,  =LCap  ci  ^  ;^LCap 

2tj2  ■'■''^12  ^12— '^12 


(Till 


Cap 


(rti 


Cap' 

123 


LCap-  _  —LCap-  ^ 

M3  '^*'13  - 

LCap-  _  —LCap-  .  =LCap-  ^L  ^  TjfLCap- 

M3  '^‘■13  '^‘■13  ^\-i  -  '^‘■13 


(rtS''' ) 


>0 


231 

Cap" 


(7232 


(  —LCap 
V/233 


(Vn) 

TT^l  <^if 

(v^) 

^23  ^^2^ 

(vi) 

^  -K2 

(v^) 

<  > 

(Vn) 

IV 

(V23) 

^^LCap^ 

—LCap^  ^  —LCap^ 

^13  ^^13 

^^LCap^ 

—LCap^  ^  :=LCap'^ 

/4.23  —  /t'23 
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^LCap  ^  _=LCap 
'^‘■12  -  '^‘■12 

—LCap^  ^  _=LCap 

^‘■13  -  '^‘■13 

—LCap-  ^  _—LCap 

'^23  —  '^23 


('74“"') 
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